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1.0 INTRODUCTION AND SUMMARY OF RESULTS 


This report presents summaries of the analytical efforts and system stud- 
ies conducted under Contract No. NAS6-2621, in connection with the GE0S-3 
radar altimeter program. The main investigative areas reported herein are in 
sea state estimation and data processing related topics. 

t 

Section 2.0 discusses significant waveheight (H-^) estimation and the 
deconvolution of the waveheight probability density function (PDF) . The al- 
gorithm presented herein differs from other estimators known to the writer in 
that the waveheight probability density function is not constrained to be a 
Gaussian. The algorithm shows good agreement with other algorithms, how- 

ever, the ensuing density functions are at times found to be highly non-Gaus- 
sian. This result applies both to use of .the standard GEOS-3 waveform data 
and use of time-realigned waveform data. The departures from Gaussanity are 
generally of the bimodal nature and, as such, cannot be adequately charact- 


erized by terms such as skewness or kurtosis. The bimodality resembles the 
theoretical solution given by Longuet-Higgins for a long-crested sea. 

Section 3.0 examines the general problem of significant waveheight esti- 
mation and presents an analysis of the resolution available from any (unbias- 
ed) risetime-based estimator. This analysis shows the H algorithm pre- 
sented here and the one currently in use at Wallops Flight Center (WFC) pro- 
vide performance which is close to the theoretical resolution limit. The 
principal inference of this result is that a sophisticated (waveform based) 
Hf/s estimator does not exist which will yield significant improvement in res- 
olution over the algorithm presently in use at WFC (developed by Hayne) . This 
analysis also quantifies the increase in performance achievable with the WFC 
estimator when it is coupled with the time-realignment technique developed 
by Walsh. It should be emphasized that the performance analysis given in 



Section 3.0 pertains to estimators based on changes in the leading edge of . 
the waveform as a function of sea state. The performance analysis work 

given in Section 3.0 also leads to the deveTopment'>of a totally different, 
low waveheight, estimator termed the variance-based H^y^ algorithm. This 
estimator is felt to offer potentially higher resolution in the- low waveheight 
range compared to other waveform estimators since use of the variance behavior 
of the ^ estimates introduces .an additional modeling element into the prob- 
lem. This estimator is presently being compared with buoy data - the results 
will be given in a later report. Appendix -B contains information on another 
Hf/s estimator which was investigated. In it, the auto covariance -of the wave- 
form plateau region was used as a wave-height sensitive parameter. The re- 
sults were negative in that the covariance was found to be only weakly sensi— \ 
tive to sea-state changes. 

Section 4.0 presents a discussion of GEOS-3 backscatter data for periods 
in which the radar cross-section (<?°) appears to 'increase markedly. These 
periods were observed early in the GEOS-3 program and were initially thought . 
to be due to anomalous scattering conditions. Comparisons of -theoretical . and 
measured (7° values indicate that these periods represent relatively ’calm or 
swell-dominated sea conditions. This section also discusses the data periods 
in -which increases in the attitude/specular gate have been observed. These . 
increased values are shown to be due to ocean surface inhomogeneities. Com- 
parisons between published 0° and those derived from the two available AGC 
calibrations are given. These comparisons show the "clutter" calibrations to 

be in considerably better agreement with published and theoretical values, 

* 

Section 5.0 covers work relating to system- and data processing considera- 
tions. The items discussed comprise altitude data editing, sea state .altitude 
bias effects, precipitation sensitivity of .the radar altimeter,- waveform sam- 
pler corrections, and tracking jitter correlation properties. 
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2.0 ESTIMATION OF SIGNIFICANT WAVEHEIGHT AND DECONVOLUTION OF THE WAVE- 
HEIGHT DENSITY FUNCTION 

• The GEOS-3 average backscattered waveform can be modeled as a repeated 
convolution of three functions, two of which describe the ocean surface fea- 
tures with the third representing altimeter effects. In mathematical form, 
the preceding may be written as follows: 

y(t) = h(t)*p(t)*f (t) (1) 

where 

y(t) represents the GEOS-3 average return waveform, 
h(t) stands for the ocean surface roughness function, 
p(t) is a function which incorporates altimeter wave- 
shape and tracking jitter effects on the return 
waveform, here called point target response, 
f(t) represents the ocean flat surface function, 
and * denotes convolution. 

Combining p(t)*f(t) into the single function u(t) leads to the expression 

y(t) = h(t)*u(t) (2) 

Given y(t) and a model for u(t), the problem is to estimate h(t), the sea sur- 
face roughness function, by performing the deconvolution specified by equation 
(2). Determination of h(t) is encumbered by the fact that y(t) is available 
only in terms of sixteen noise-perturbed samples with nominal 6.25 nano sec. 
spacing and because u(t) is not precisely known. Given the surface roughness 
function h(t) it is possible to deduce significant waveheight, H.^^* by inter- 
preting h(t)‘ as a surface roughness probability density function. 

An attempt has been made to formulate and solve this problem in a general 
manner. To this end the following developmental guide lines were followed: 
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(a) Points in time at which y(t) is sampled are no'f constrained 

(b) The only restriction placed upon h(t) is' that ‘its Laplace trans- 
form, H(s), be rational in s. ■ ‘ •' • '• “• " v " 

In the sequel the following steps will’ be taken so- that, using the model de- 
fined by equation (2), h(t) may be determined from experimental- data: 

(a) Apply the Laplace transformation to equation (2) and linearize 
the resulting expression using quasilinearization. 

(b) Inverse transform the linearized model and discretize the resulting 
time-domain equation in order to account for the sampled nature of 
the problem. 

(c) Develop equations which enable specification of the constants in 
the rational function H(s) from measurements of y(t. ). 

• r iC 

(d) Utilize h(t) to compute significant waveheight. 

2.1 Linear System Model and Algorithm for Determination of h(t) 

As already mentioned it is necessary to approximate u(t) which 'is used 
to model the effects of the GEOS-3 altimeter point target response and the 
ocean flat surface function. In the present study, the ocean flat surface 
function is represented by a unit step function with transition at the time 
origin. On the other hand, the GEOS-3 altimeter point target response, which 
has, roughly speaking, a Gaussian shape, is approximated with a truncated, 
raised cosine function. Thus 


U(s).= L[u(t) ] 



X 
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T = one half the pulse width of the altimeter point 


1 

s 


target function 

represents the transform of the ocean flat surface 


response function, and 


2 . 2 
s +w 


is the transform of the point target function approximation. The general 
shape of u(t) is illustrated in Figure 2.L. In equation (3), T is selected 
so as to reflect the entire altimeter point target response function, includ- 
ing jitter, according to the expression 

O — 0.362 T 
Pt 

where cr is the standard deviation of the altimeter point target function, 
pt 

Taking. the transform of equation (2), the time-continuous GE0S-3 return 
waveform can be expressed as 

Y(s) = H(s) U(s) 

= (N(s)/D(s)) U(s) (4) 

where N(s) = a s n + a n s n 1 +. . .+ a 

o 1 n 

D (s) = s 11 f b 1 s 11 " 1 +. . .+ l n . 

The problem is to determine {a_.,K} . In order to obtain a computational al- 
gorithm for solution of this problem, it is convenient to convert the highly 
nonlinear form represented by equation (4) into an iterative linear problem. 

Applying quasilinearization [1], a Taylor-series-like functional expan- 
sion, to equation (4) results in the following relation 
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Figure 2.1. Illustration of general shape of 
system model input function. 
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( 5 ) 


Y ±+1 ( S ) = Y i (s) + (N 1+1 (s)/D i (s)) U(s) - (D i+1 (s)/D i (s)) Y*(s) 
mm 111 


where 


..1 , . i n-1 i n-2 , i 

N (s) = a^s + a 2 s + — + a^ 


D^Cs) = s n + +. ..+ b* 


Y 1 (s) = model of GEOS-3 average return waveform 
m 

tti 

determined for i — iteration 


i = iteration number = 0,1,2,.. 


As shown in reference [2], the model response function Y^ (s) can be express- 
ed in the time domain, using vector-matrix notation. For example, applying 
the referenced procedure to equation (4) leads to the following result 


where 


x(t) = Ax(t) + Bu(t) 
y(t) = PpX(t) + a Q u(t) 


A = 


-b 


n 


n-1 


-b 


n-1 


~b. 


( 6 ) 


B - (0,0,. 


P — j (a — b a), (a - b .. a ),..., (a- b- a ) 

p L n n o n-1 n-1 o 1 1 o j 


I , = (n-1) x (n-1) identity matrix 
n-1 

x(t) = time derivative of system state vector, 
and prime denotes matrix transposition. 

In developing equation (6) from equation (4) , initial conditions are assumed 
to be zero. Applying to equation (5) the above procedure which led to 
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equation (6) results in the following representation of model output: 


= *i<*> >r + x ii (t > p r - p b 


,i+l 


i+1 


( 7 ) 


i+1 

with P = (a 


i+1 i+1 


* Vr 


aj + V 


i+1 _ i+1 i+1 x+1 , 

P b " ( b n » Vl* ’ b l } 


i+1 x 

x^(t) = state vector corresponding to H^(s) = N (s)/D (s) 

i+1 i 

Xj^Ct) = state vector corresponding to H^Cs) = D (s)/D (s) 


Denoting the average observed (i.e. experimental) GEOS-3 return waveform 
by y Q (t), the error between model and observation can. be written as 


i+1 


i+1 


e (t) = y (t) - y (t) 
o m 


( 8 ) 


which can be interpreted as shown in Figure 2.2. From the figure it is noted 
that Hj.Cs), corresponding to Xj.(t), is forced by U(s) while Hj.j.(s), correspond- 
ing to Xj-j-Ct), is forced by Y^(s) . Selecting an index of performance for 

i+1 

objectively evaluating the goodness of fit provided by y (t) to y Q (t) and 
solution of the resulting optimization problem enables one to determine the 
{a_. , b_. } . In obtaining the solution note that at convergence (i.e. | £ -£ |-^0) 

equation (7) reduces to 

y (t) = xl(t) P 
m la 

Thus at convergence Hj.(s) represents the desired system function while the 
effect of Hj.j.(s) vanishes from the model. 

Since y (t) is known only for t = t (k=l,2, . . .,16) , discretization of 
equations (7) and (8) is required. Thus 
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(9) 


p r + p r - <ii<*k> p i 


,i+l 


,i+l _ i 


£±+1( V - W - 4 + X • 


( 10 ) 


k — 1,2 , ... ,16 


Assuming that u(t) and y m (t). change linearly between t^ +1 and. t^, it. can 
be shown that x^.(t^) and (t^) in equation (9) can be written as follows: 

VW “ ® i( V x i (t k ) + a X E »V 


+ (A 1 ) 1 Ee i <A k ) - A^l] E 


x II (t kH 5 $ ( \ )x II ( V + 0 C V B y m^ t k ) 


+ (A 1 ) 1 [0 ± (A k ) - A^.1 j B -i- 7 g (tk+1 ^ 


where 


\ t k+l ” fc k 


$ (A^.) - state transition matrix corresponding to 

, i+1 ± 

H z (s) = N (s)/D 1 (s) 


(ID 


( 12 ) 


0 i (A k ) = (A 1 ) 1 [$ i (A k ) - Iq ] 


0 


-b J 


n 


n-1 


■h 2 - , i 

b n-l **• " b l 
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and B is as previously defined. Note that the assumption of linear variation 


of y m (t) between and , while it might affect convergence rate will 

not influence the approximation of y Q (t). ~ 

Combining equations (9) and' (10) results in 'the vector error expression 


e 1 " 1 " 1 = Y - X P i+1 
• o — 


(13) 


where 


X = 






x ll(tl) 




Y = 


W + x i’i (t i ) p b 


7 o (t 16 ) + x ld]d p b 


.i+1 


r p i+i 

_a_ 
„i+l 


i+1 

and e is a column vector of errors between observation and 
model . 

Selecting the performance index as minimization of the sum squared error , 

. i+1 * i+1 . , . , . . 

(£ ) • e , results m the parameter estimate 


where P 


i+1 


pi+1 


/\. 

p i+l 

a 

pi+1 

^b J 


= (X*X) ■ L X’ Y 
' — 0 


(14) 


Given initial values for {b^ } , it is feasible to evaluate equations (9) and 
(10), using (11) and (12), and to update the starting estimate of the unknown 


parameter vector, P , by evaluating equation (14) . 
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The computational algorithm for evaluating equation (14) can be initi- 
alized by specifying {b?} and y°(t. ) or, equivalently, by specifying {a?,b?}. 

1 m k. - . 4 4 1 

The initialization procedure that was used to compute the results presented 

in paragraph 2.4 consists of setting {b°}. equal to coefficients which define 

3 - ' 


the all pole maximally flat, unit delay approximation discussed in [3]. Also 
since y i+1 (t,) must approximate y (t, ), the iterations are started under the 

HI K O K 

conditions y (t, ) - y (t, ) until the computation stabilizes at which time use 
m k o k 

i 

of y m (t^) can be made for the (i+l)th iteration. Stopping criteria for the 
algorithm described above are readily formulated. A simple stopping criteria 
would be 


(e i+1 )’ e 


i+1 


- (eVe 1 


> 6 , continue 
< 6 , stop iterations 


where 6 is an appropriately chosen small positive number. 

2. 2 Alternate Interpretation of Model 

By manipulating equation (3) the model developed in' paragraph 2.1 can 
be modified such that ideal altimeters (i.e. point target response equal to 
the Dirac delta function) may be represented'. Consider the limit of equation 
(3) as T 0 , i.e. 

U (s) = lim U(s) 

T>0 

“7 * (15 > 

Therefore, all of the equations developed in paragraph 2.1 apply to the pres- 
ent' situation provided 

u(t) = Vl (t) , 

that is, u(t) is set equal to the unit step function u(t). This form of the model 

' * ’ * " ' 2 

has been used and results from its application will he presented in paragraph 2.4. 
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Another interpretation results upon combining equations (3) and (4) to 


write 

Assuming zero initial conditions 

^ - i l [«w] 



By normalizing dy{t)/dt such that the curve defines unit area, this result 
can he interpreted as the surface roughness probability density function (it 
is assumed that y(t) is monotone nondecreasing). A somewhat similar result 
follows if the altimeter point target response is approximated with a rec- 
tangular pulse. 


2.3 Analytic Evaluation of Surface Roughness from H(s) 

The function h(t) (or h(t^)) can be interpreted as the unnormalized sur- 
face roughness density function. Define 


h(t) = f h(t) 
a 




where 


a - 


h(t)dt 


and Tg is the time expanse of the GEOS-3 waveform sampling gates. 


Then 


a 


2 



(16) 
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where 



T g 

fck h(t)dt 
o 


In the context of paragraph 1.1 


P' A -1 
a 

m 2 a 


{ 


.-1 


$ (Tg)Tg - 2TgA (Tg) 
+ 2A _1 G(Tg)J. B 


and 


"l = 


P* A 
a 


-1 


\ $ (Tg)Tg - 0(Tg) 



2.4 Computer Implementation and Simulated Results 

The algorithm for determination of H(s), presented in paragraph 2.1, has 
been programmed and tested using GEOS-3 frame-averaged gate data. Details 
relating to use of the computer program are contained in Appendix A . Appli- 
cation to computation of surface* roughness density, h(t), is the topic con- . 
sidered in paragraph- 2 . 4 . 1 while. its use as an estimator' of is. treated 

in paragraph 2.4.2. ' . - . 


2.4.1 Computation of Surface Roughness Density Function 

A 

The algorithm developed in paragraph 2.1 computes h(t), the surface 
roughness density function, given the frame-averaged GEOS sample gate data. 

This is an inverse problem for which numerical computations can be unstable 
due to noise and modeling error effects in the solution technique employed. 

The method used herein avoids explicit numerical deconvolution since h(kT), 
k = 0,1,2, ... ,N, is not computed directly in terms of the measured or experi- 
mental data. Rather, h(t) is obtained (in terms of a small number of parameters) 

14 



by applying system identification techniques. 

Typical results obtained from applying the algorithm will be presented 
at this time. All- results shown were obtained with u(t) (see equation (2)) 
defined as a linear ramp starting from zero at gate one, extending over four 
gate intervals, and equal to unity thereafter. Figure 2.3 shows the sampling 
gate waveforms for three consecutive frames of Rev. 6893. In this figure,, 
the vertical scale reflects the measured values of frame 103 and frames 104 
and 105 were adjusted such that their normalized response asymptotes (see 
discussion in Appendix A ) were equal to that of frame 103. This enables a 
relative comparison of the rates at which the three response curves rise. 

The deconvolved surface roughness probability density functions corresponding 
to the return waveforms of Figure 2.3 are shown in Figure 2.4. Note that the 
tendency of frames 104 and 105 to rise early relative to frame 103 is reflect- 
ed in the probability density curves of Figure 2.4. In particular the two- 
step shape of frame 105 (Figure 2.3) results in a bimodal density function. 

The ^ values shown in Figure 2.4 were computed using Hayne’s [4] algorithm.' 

values computed from the probability density curves shown' are placed in- 
side parentheses. 

Bimodal densities have been observed in a number of instances for large 
sea state conditions (H^y^>7.0 meters). It is natural to question whether 
or not this shape might be due, at least partially, to tracking loop jitter 
effects. Figures 2.5 and 2.6 present results obtained from processing both 
time realigned (dashed curves) and unrealigned data (solid curves). In the 
case of Figure 2.5, time realignment resulted in the attenuation of the ten- 
dency toward bimodal behavior. In contrast, however. Figure 2.6 shows that an 
approximately trapezoidal shape reverts to bimodal nature upon application of 
time realignment. This result is interesting if it is noted that after time 
realignment both the shape of the density function and the value of are 
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NORMALIZED RESPONSE 



SAMPLING GATE NUMBER 

Figure 2.3. GEOS-3 Sample Gate waveforms corresponding to the surface 
roughness density functions illustrated in Figure 4. 



REV 


6893 


FRAME 

103 

H 

1/3 = 8.3 (8.1) 

FRAME 

104 

H 

1/3= 7.5 _ (7.5) 

FRAME 

105 

H 

1/3 = 8.1 (8.2) 



to the waveforms of Figure 4. 



REV 6893 , FRAME 107 



SAMPLING GATE NUMBER 


Figure 2.5. An illustration of the effect of time realignment on 
the surface roughness density function. 




SAMPLING GATE NUMBER 


Figure 2.6. An illustration of the effect of time realignment 
on the surface roughness density function. 


more consistent with results from neighboring frame data than is the unre- 
aligned data. 

A 

In its present form the algorithm cannot cope with h(t) functions en- 

A 

countered when H^^<4.0 meters. In this range of , ‘sea- state values h(t) be- 
gins to approach the Dirac delta function (in the limit as it is 

equal to the delta function) which cannot be approximated by a finite rational 
function. Figure 2.7 is an illustration of this effect. Note the increased 
tendency to oscillation displayed by the density function plotted. This be- 
havior is a result of the approximation problem mentioned above. 

Another characteristic of the approximation technique is that the shape 
of h(t) is dependent upon the order chosen for the approximating system. For 
convenience let the respective orders of numerator and denominator polynomials 
of H(s) be denoted by (m-l,n) where 

0 < m-1 < n 

1 c n £ 5 

Figure 2.8 illustrates the effect of model order upon the resulting probability 
function. In this figure, the notation (3,4), (2,3) should be interpreted to 
mean that the resultant curve was obtained as the arithmetic mean of two iden- 
tification procedures - one for which H(s) had polynomial orders (3,4) and 
the other for which the orders were (2,3). The example in Figure 2.8 shows 
that the average density produced by (3,4), (2,3) is somewhat different from 
that realized when (2,3), (1,2) is used. All of the curves presented in this 
section were computed as an average of two model fits to the sampled gates. 

In all cases, either (3,4), (2,3) or (2,3), (1,2) combinations were used in 

i 

A 

obtaining h(t). Use of the combination (m-l,n) = (4,5) frequently resulted in 

A 

highly oscillatory h(t) and experience indicates that it can only rarely be 
used. 
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Figure 2.7. Estimated density function for H. 
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SAMPLING GATE NUMBER 


Figure 2.8. Illustration of the effect of model order upon approximation 
to surface roughness density function. 


2.4.2 Estimation of Sea State 

The algorithm developed in paragraph 2.1 can, of course, be used to esti- 
mate sea state for all values of h -j _/3 if u(t) in equation (2) is made equal 

, t * 

to the unit step function. This redefinition of u(t) avoids the approxima- 
tion problem discussed in paragraph 2.4.1 as K^-J-O. Therefore, the results 
presented in this paragraph were obtained under the following assumptions: 

(a) u(t) is the unit step function 

(b) Only frame averaged (3.2 second average) data 
is to be processed 

(c) (m-l,n) = (2,3) was used in obtaining all 
results presented. 

Given that u(t) is the unit step function in equation (2), may be comput- 

ed from the expression 


Kl /3 = °* 6 Vl ff2 -°p t l * Signum(a 2 -0^)- 


where 


and 



is given by equation (16) 

= altimeter point target function, 
(including jitter) 


(17) 


A sample of the computational results obtained from simulations is pre- 
sented in Figure 2.9 where they have been compared, in the form of a scatter 
plot, with the results produced by Hayne’s algorithm [4]. It is noted from 
the figure that for H meters, the two approaches produce results that 
are generally in good agreement; however for smaller the disparity is 

significant. 


2.5 Development of a Constrained Estimator of. Significant Waveheight 

This paragraph describes a method which provides least square estimates 
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Figure 2.9. Scatter plot comparing present algorithm with that of [4]. 
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of under the constraint that H^y ^ 1 0. Development of the algorithm is 

made under the assumption of Gaussian statistics. Utilizing notation intro- 
duced in paragraph 2.1, the model of the GEOS-3 average return waveform can 
be expressed in transform notation as 


Y (s) = H (s)U(s) 
m n 


(18) 


where 


- ST'i 


, n 


V s) - 


1 - e 


- sT 


ST 


by definition. 


T q = a time-domain shift parameter, 

X = parameter related to surface roughness.- 
th 

Here ^(s) represents an n= order convolutional approximation to a Gaussian 
probability density function as discussed in [6] with n=l,2,... . The basic 

function used in generating the Gaussian density is a uniform density of width 
T nano sec. with leading edge displaced T q from the origin. For U(s) as spec- 
ified in equation (3) of paragraph 2.1 and with n = 3, the discretized model 
output in the time domain is given by 




i=0 

3 

-2 

i=l 


U(z ± ) 
2T T 3 


r 4 


r 3 

- 

'l 

ft 

24 

1 ’ 
2 

a + 

— — (co z . - Sin z . ) 
3 x x 


L 

a) 

— 

to J 


U(w ± ) 

3~ 

2T T 


r 4 

r 3 -i 


) W i 1 

) 24 2 

L w 

w i l ' 

— : x- (tow . - Sin to w . ) 

6 3 x l 

L. CO -j 

] 


(19) 


where 


z . = 


w. 

x 


t - T - i T 
o 

t - 2T - T - iT 
o 


, t = 0,1,2, ... ,15 


p(*) = unit step function 


0 ) = — 
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and T is selected to approximate the ■ altimeter .point, -target response function. 
Note that in the discretization specified^ in equation (19-) that the time axis 
is scaled in units of 6.25 nano sec. 

Defining the approximation error as 


2 

e = 


2 


( y o (t) ‘ y m (t) ) 


( 20 ) 


5% 2 a 

the problem is to find T q and T which minimize £ such that T ■> 0. 

For this particular formulation it is straightforward to find approxi- 
mate T q * and T* by implementing an exhaustive two-parameter search procedure. 
Given t*, the standard deviation of surface roughness can be calculated as 


a = 6.25 (~) 


nano sec. 


Therefore, significant waveheight is given by 


H l/3 * 6 ° 


= 1.875 T* meters 


A computer program-based upon the procedure described above was written 
and used to compute from frame averages of GEOS-3 return waveforms. Be- 
fore computing x*, the return waveform was preprocessed as follows: 


y a - (i V «) /4 

v t=0 ' 


y o = 


16 \ 

E y 0 (t>) 

t=14 7 


y o (t) * [ y 0 (t) ■ y ]/[ y o- y o] 


/ 3 


t = 0,1, . . .15 
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where <■ is used to indicate a redefinition of y Q (t) . This preprocessing is 
performed so that the fitted function approximates a cumulative distribution, 
but could.be avoided by introducing an unknown multiplier into equation (18) 
and performing a three-parameter search. 

Computational results from application of this algorithm are presented in 
Figure 2.10 which again makes use of results obtained from applying the tech- 
nique described in [4]. So long as H^y^, as computed using the method of [4], 
is greater than about 2 meters the two results are in general agreement. Again, 
for < 2 meters the correspondence is not good. 

2.6 Summary and Conclusions' 

The results obtained from application of the algorithms developed herein 
lead to the conclusion that linear system theory concepts can be successfully 
applied to significant waveheight est ima tion. If the estimated sea state ex- 
ceeds approximately 2.5 meters, both of the algorithms presented herein pro- 
duce results that are in satisfactory agreement with the method developed by 
Hayne [4]; however, for less than 2.5 meters erratic performance of the 

estimators was noted. It will be shown in Section 3.0 that estimation of sea 
state for calm sea conditions (i.e. Ry^ ^ 0) is a very unstable problem for 
which the variance of any unbiased estimator can be very large, for data 
rates prescribed by GEOS-3 operational parameters. 

Attempts to determine surface roughness probability density functions 
has been successful provided that the associated H.y^ > 4.0 m. Since the 
algorithm can be extended so as to cope with instability encountered as 
H l/3 ttie a PP roac h used appears to be capable of describing surface rough- 

ness probability density for a wide range of H_y^ values. 

The algorithm of paragraph 2.1 has been used to compute an -d com- 

parison shows good agreement with results obtained from an existing method 
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Figure 2.10. Scatter plot comparing constrained 
estimator with that of [A]. 
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[4] so long as > 2.5 m.~ 

The constrained estimator of paragraph 2.5 performs satisfactorily and 
is intuitively appealing since negative H^y^ is physically meaningless. How- 
ever, the important question of bias induced by this- estimator is a matter of 
concern. By resorting to Monte Carlo simulations, this question could be re- 
solved but such a study falls outside the scope of the present investigation. 
This algorithm might he generalized to enable fine structure identification 
by utilizing the sum of Gaussian-shaped densities . In this more general set- 
ting it is doubtful if a multidimensional search approach as employed here 
would remain a viable computational approach. : 
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3.0 ASSESSMENT OF PERFORMANCE LIMITATIONS OF GEOS-3 SIGNIFICANT -WAVEHEIGHT 

ESTIMATORS 

A number of algorithms have been developed ' for estimating ocean signifi- 
cant waveheight by analysis of the GEOS-3 altimeter average return waveform. 
The results of limited comparisons of three’ such algorithms is presented in - 
Section 2.0 of this report, from which it is concluded that acceptable agree- 
ment exists between the algorithms for greater than about two meters. 

However, if the estimated is less than two meters the estimators all 

show that H.^ resolution degrades rapidly as -> 0 (this result was antici- 

pated by Miller and Brown [1]). The purpose of this report is to evaluate the 
performance limitations encountered by any estimator of which uses the 

GEOS-3 models for average return waveform and noise. The Cramer-Rao inequali- 
ty will be used to establish estimator performance bounds. 

The objective of this study is to -determine estimator performance 

limitations and to qualitatively assess .the degree to which the present GEOS-3 
H^ 3 algorithm [2] approaches the resulting bound. The Cramer-Rao bound is 
selected for this investigation because (1) compared with other bounds (e.g. 
Baranken, Ziv-Zakai, etc.) it is easy to apply and (2) the resulting bound is 
an upper bound on performance in that other tighter bounds show that this up- 
per bound cannot be achieved in practice (see Seidman [3])- In the above con- 
text, then, the Cramer-Rao bound represents a relatively easily applied, yet 
severe, test of algorithm performance. 

3.1 The Cramer-Rao Inequality 

Development of the Cramer-Rao inequality is readily available (see for 
example [4], [5] and [6]) and states that any unbiased estimator must satisfy 
the following relation 
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[(e-9) 2 ] > 


[jq fnpCx^x^ . 


.x 


i e i> 6 */ 


J 


( 21 ) 


where 02 =s ®2^^ * a •^ unc ^-*- on °f ® 

0 and 0^ are unknown constants to be estimated from 
experimental data 

p (x^,x 2 , . . .x^ [ 0^ 9 8^) denotes the probability density 
function of measured data x^x^.-.x^^ given 
9^ and 0^ 

✓ N . 

8 is an estimate of parameter 0 
and E denotes expectation. 


In this analysis p(x-,x , ...x [0,0) is assumed to follow the Gaussian law 

±2 n' 1 L 

since the experimental data from which significant waveheight is estimated is 
obtained by linearly combining a large number of individual noisy return sig- 
nals. 


3.2 Analysis 

For purposes of this analysis the ideal 
turn waveform is assumed to be specified by 

y(t) = y(t,e s e ) = - 

/2? 


, normalized GEOS-3 average re- 
[7]* 
v 



( 22 ) 


—00 


where v = (t - , 

0^ = constant , 

9 2 ' °tt + < H 1/ 3 / 2c > 2 ■ 

2 

represents altimeter point target effects , 


*A similar analysis is given in [7], however, 
to rise-time not was used. 


the Cramer -Rao bound relating 
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is significant waveheight , 

and c is' r the- speed of - light in units of meter /nanosec. 

The time-averaged GEOS-3 return waveform at the output. of the sampling gates 
is represented by the following model [7] 

x(t ) = y(t ,0-,0-) + z(t ) , n = 1,2, . . . ,N (23) 

n n x z n x 


with 


cr?(t ) = Var[z(t )] = -4- 

z n n 




S/N = altimeter IF signal to noise ratio 
2 

F = a system constant which for a 1-second averaging 
period is equal to 200 
N = 16 

x(t^) = observed signal, at n, th sample gate 
yCt^, 0 ^, 62 ) = observed signal in absence of noise 
z(t^) = additive, independent Gaussian noise 
t -t = 6.25 nano sec, n = 1,2,... EL 

In this model y(*) represents the true received waveshape and z(*,) the noise 

which arises mainly from the fluctuating nature of the received signal. For 

substantial averaging periods z(*) is Gaussian by the central limit theorem. 

The present study will be concerned only with the high signal-to-noise case 
_2 

so the term (S/N) in the noise model will he dropped. 

For the above model assumptions, the likelihood probability function for 
observing x(t^) ,x(t 2 ) > ♦ • .x(t^) given 0 ^ and is, since the noise is inde- 

pendent and Gaussian distributed 
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N„ 


L (x^ , x 2 , • . • x^ [ *^1/3^ 


K '1T 5"ul ^ 

n=l z n 


f- [*(t n ) - y(t n )J 


- 2 


2 At ) 

z n - 


Therefore, the likelihood equation is 


N i N i r" i 2 

V -v-i [x(t ) " y (t ) 1 

1 <V‘2’”' 3 ‘J 0 l’ H l/3 ) -* 1 - 2 . - 2 ^73 ■ (24) 


n=l 


n=l 


z n 


where = £nK = constant. 


Taking the partial derivative of 2, with respect to model parameter 

(q i = e i ’ q 2 ‘ hrf sives 


a a 

8q 


N i 

_ V 1 ag z (t n ) 
1 " 2 1 W 3q i 


+ i y pK>- y(t n )] ay(t n ) 

i=l L 


oht) 

z n 


3 q. 


Tx(t ) - y(t jT a At ) ^ 

i_ n n J z n 

3q 


cA t ) 

z n 


> ; i = l,2 . 


(25) 


After taking the derivatives indicated in equation (25) , substituting for 

2 2 4 

a (t ), squaring the result, neglecting all terms divided by F or F and 

Z XI 

taking the expectation there results 


(M.\ 

“i 

V 1 I 

f »y(t n ) \ 

_ 

' 2 a 2 (t ) 1 

n=l *7 

V 3q i / 


i = 1,2 


(26) 


From section 3.1 the Cramer-Rao inequality states that 
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(27) 


r. N. 


E '[ < V«i )2 ] 


IT 


V ' 1 

^ a 2 (t ) 
n^l z n 


'3y (t p ) 
3q i 


2 “I" 1 


J",CV.' SSi": \ 


where q^ = 0 and q 2 =H l/3 

. . 

r\ J" ,* , * 

Equation" (2’7) can'.*be“optimistic when more than one' parameter must be estimated. 

< ~ •>* „ I * 

For the case at hand both 0^ and “wst be determined; therefore, as dis- 

. * .v 

cussed in [6] the Fisher’s information matrix is appropriate and leads to the 
following definitions 


J = ( J . . ) = Fisher Information matrix 
ij 


with elements 

; i» j - 1,2 . 

In the multidimensional case,Tthe Cramer-Rao inequality takes the following 
form [6]: 


J. . = E 
ij 






Sq. 3qj 


E K-‘ , i )2 }'- ( j ’ 1 )ii ; i -* i:;2 - • (28) 

With the aid of Leibnetz’s rule the derivatives required for the Fisher in- 

V * . . 1 * . « ' f . 

formation matrix can be evaluated and are 



100 XT’ 

F8 2 Z, 

n=l 


-.w , 


y(t n ,e r 0 2 ) 


(29a) 
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1/3 


1/3 j 
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J 12 -e 


( s“ ) ( 4*3 )] 


N. 


- 25 H 


= J 21 = 


1/3 


Zend; 


1 - ■- . • • -w 

<W * 

y«n, e r e 2> 


(29c) 


n=l 


where w = 


[ 2c(t n- e i^] 

(2c0 tt ) 2 + *y 3 


For J ^2 ^ 0 (i.e. non-zero correlation between 0-^ and-H.^^) It follows that 

A 

(see [6]) the variance of the estimator h ^/3 satisfies the expression 


(X* 

°H 


tt(2c) 4 02 


1/3 


22 100 H 


1/3 


*1 

2 

n=l 


^WV 


- w 


(30) 


As note the asymptotically unbounded nature of estimator variance 

Cf£ 2 . Thus for calm sea conditions, any minimum variance unbiased esti-? 

“l/3 

mator can be expected to display poor performance. The three estimators dis- 
cussed in Section 2 of this report are characterized by erratic performance 
as H^y^-J-0. In Section 3.4 the Cramer-Rao bound for GEOS-3 H.^ estimators 
will be presented from results obtained via a numerical evaluation of the 
inverse of Fisher’s information matrix, equation (29). 


3.3 Modeling Altimeter Tracking Loop Jitter 

In the foregoing analysis the model employed in the development leading 
to equations (29) did not consider altimeter tracking loop jitter effects. 

Since this is an important effect the resulting Cramer-Rao bound must reflect 
its presence. Brown [8] studied jitter effects on the Skylab S-193 altimeter 
performance. In this section results from the analysis in [8] will be adapted, 
using approximations, for use in the present study. In effecting this 
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adaptation the goal is to achieve a closed form, easily evaluated result from 

* 

inductive arguments based upon assumptions. ' 

Tracking loop jitter manifests itself in two distinct ways: 

(a) Introduction of a smearing effect on the average ideal 
return waveform y(t, 0 ^,& 2 )» an -d 

2 

(b) Enhancement of the altimeter noise process a (t) 

z 

With r-especb to (a), jitter will be modelled as an effectively increased 
point target .function which combines in a root sum square sense.' - That is, 
with jitter included, the point target function is 


°tt). * («tt + ) 


1/2 


nano sec. 


(31) 


where = 4.0 as determined from experiment. 

Noise enhancement,' (b) above, as shown by Brown [83, can -be- described 
by the .following convolutional sums (using notation of this chapter-)--* 


-E[z 


(t)] * K o 2 P m y(t +mT) 


(32a) 


m =_oo 


and- • 


E[/(t)] - K* ^ P m {[y(t+mt)+ (|)] 

m=-®° * - 

+ y^ft+mx)^ .*(32b) , 

where " K = "a system constant 

E is expectation operator 
P m “ probability masses of the discrete 

, 4 r 

jitter probability density function 
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S/N = IF signal to noise ratio 

* 

T " time quantization, parameter of the 
altimeter tracking loop % 
z (• ) = average return waveshape with 
jitter included. 


Assuming for the moment that a hypothetical radar average return pulse is a 
step function (i.e. transmitted pulse is the Dirac delta function), that time 
quantization is small (i.e. T->"0), and S/N->-» equation (32) can be expressed 
as follows 


E[z(t) ] = K q 


CO 


I 

~oo 


p m (v) y (tW) dv 


(33a) 


E[z^(t)] 


(X 

■<J 


•2p (v) 
m 


y^(t+v) dv 


(33b) 


with y (•) a unit step function. 

Thus the average noise effect of jitter, for the special case considered, can 
also be described by a spreading effect on the average return waveform coupled 
with a nonlinear combination. Assuming an identical phenomenological noise 
effect in the non-ideal GEOS-3 altimeter results in an increased noise level 
given by the relation 


a z(t) = E t z2 < £ )l - [Ez(t) ] 2 

= 2y(t,0 1 ,6 2 ) - y 2 (t,0 1 ,0 z ) > 0 

where now © 2 = + (h 1/3/ / 2c ) 2 • 


( 34 ) 
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3 . 4 Computational Results 

To evaluate the Cramer-Rao bound a computer program was written and used 
to compute the elements of Fisher’s information matrix (equation (29)). The 
results are shown in graphic form in Figure 3.1 where three curves have been 
plotted with (1) jitter effects totally omitted, (2) partial jitter effects 
resulting from effective increase in point target response, and (3) including 
total jitter effects (i.e. both effects (a) and (b) of Section 3.3). It is 
emphasized that the est im ator performance curves plotted in Figure 3.1, under 

5, '• {’ ■ . ■ ‘ - . . i. " * i ,.S~- ■ 

the model employed, define performance limits which cannot be surpassed by any 
unbiased estimator. Also, the severe degradation of estimator performance for 
H^3<3 meters is significant and has been observed in the analysis of GEOS-3 
data. 


In obtaining the data for use in Figure 3.1, 0^ was used as a parameter 

to verify that the resulting bounds were, for .all practical purposes, not af- 

/ 

fected by the position, in the GEOS-3 sampling gate set, of the return signal. 

Figure 3.2, shown fox comparison with Figure .3.-1, illustrates two Cramer- 
Rao .performance limits curves, for estimators since jitter effects --were; ■«“ 


ignored altogether^and it was further assumed, that^y/.t ,-0 , 8^') /was available,.. r 6 
in cojitinuo.us form.. The dashed curve-of this. figure is an estimate of. per-—-’ 
formance achieved .by.. fhe estimator presently ..used for GEOS-3 computa- - 

tions [2]. Referring to Figure 3.1, it can be seen that the GEOS-3 ,H, ' est-i-r 

mator closely approaches the Cramer-Rao bound. 


The GEOS— 3 performance limit curve shown in Figure 3.2 was drawn from the 


relation 


E [( fi l/3- H l/3) 2 ] 


64/F c 4 Q 3 J 2 


N 


H. 


1/3 


(35) 


which following [4] was derived under the assumptions of continuous measurements 
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Figure 3.1. Calculated GEOS-3 H , resolution . 
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Figure 3.2. Cal'culated and measured GE0S-3H1/3 resolution 
compared to the calculated SEA SAT resolution. 

40 



perturbed by white noise with cr^ = 0.06. SEASAT performance .limitation was 
obtained from equation (35) assuming a point target .standard deviation. of 
3.0 nano sec. and a pulse repetition frequency ten times that of GEOS-3. 


3.5 Variance-Based Algorithm 

The above discussed bounds indicate an interesting paradox - variance of 
the estimate appears to be a very sensitive parameter for values between 

0-3 meters. This behavior suggests the possibility of merging an estimate 
based on waveform data with one based on variance behavior. This possibility 
is explored in the next few paragraphs. 

Computation of waveheight data based on variance or standard deviation 
values was first suggested by the waveheight resolution analysis given in [1]. 
For this reason, the algorithm to be given is based on a curve-fit using the 
functional form 


H. 


1/3 


.042(7.66 IL^+T 2 




+ .546 


o 2 (; 

_JLl 


7.66 H^+T 2 


) 


H 1/3 R 


where 


= significant waveheight in meters 
Cl. = altitude tracking jitter in n.s. 

T = 3 dB post detection pulse width in n.s. 
t = smoothing interval in sec. 

R = variance reduction factor 

The above equation was written as H Versus variance and the R factor em- 
pirically determined by curve filtering to a scatter-plot of values (as 

determined by the waveform algorithm) versus NOAA/SMG ground-truth data. The 
numerical values of the parameters were O ^ =3 n.s. and T = 10 n.s. The fitted 
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curve was previously shown in Figure 3.2. 

Computation of waveheight values using the variance relationship involves 
the following steps: 

1. Six frame estimates of mean and standard deviation 

are first computed using per-frame numbers (in 

meters) as input values 

2 . The standard deviation (Cf) value is then tested to 

determine if it is > .815; if not the 6 frame mean 

-■ : * -v 

values computed in the first step are used as 
values 

3. If a > .815 the variance algorithm is used in the forpi 

i 


"i/s' 



O 2 -20. 


72 


)- v<? 


63.03 a -20.72 


) 2 


370.59 



A comparison of the H^y^ values obtained using the two algorithms is 
shown in Figures 3.3 and 3.4. In each figure, H^y^ values obtained using the 
waveheight algorithm is shown as the solid line and values from the variance 
algorithm is shown as the dashed line; any discontinuities in the dashed line 
indicate that the computed waveheights exceeded 3.5 meters according to the 
condition O < .815. These figures also show NOAA/SMG estimates of waveheight; 
additionally, Figure 3.3 gives the value measured by the laser profilo- 

meter. Examination of the data shown in these figures indicates the variance- 

derived estimate to be in better agreement with the available ground truth 

& 

datal - The- rapid changes in the dashed curve hear the end’ of the data-span in 
Figure '3.3 is thought -to be due to non-stationarity of input values' rather 
than to real waveheight changes. Additional ground truth data is needed to 
fully 'evaluate the variance algorithm* --note that the one available measured 
value (from the profilometer) is in very good agreement with the algorithm. 
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Figure 3.3. Comparison of Hj ^/3 values obtained from the WFC algorithm (solid curve) and from the 

variance-based algorithm (dashed curve), SMG and profilometer values are also shown. 
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Figure 3.4. Comparison of H 1/3 values obtained -from the WFC algorithm (solid 
curve) and from the variance-based algorithm (dashed curve) . 



3.6 Conclusions 


The GEOS-3 estimator currently in use closely approximates the per- 

: '< >1 ■ .< ■ . . : 
formance bound developed from the model presented, computation is a 

highly nonlinear problem which for GEOS-3 is solved by conversion to a linear- 
ized model which is iterated to convergence [2]. It has been shown that such 

; J « 

linearized estimators, under ■ appropriate conditions, are unbiased [9,10]. 

This implies that the current GEOS-3 estimator, while a suboptimal one, 

achieves near-optimal performance. 

A significant result of the analysis is that under calm-sea conditions 
(i.e. H^yg-*0) the performance of any estimator of can be expected to 

exhibit marked degradation. Two obvious techniques for combating this prob- 
lem are (a) use of higher pulse repetition frequencies and (b) reduction of 
point target effects, a . Of course neither of these options is applicable 
to GEOS— 3. 

There is a theoretically optimum H i n the sense that estimates of 
this one particular value will have minimum variance compared to that of all 
other values of For the GEOS-3 radar parameters, this optimum value 

is within the 4.0 - 5.0 meter range and is characterized by a very broad mini- 
mum (Figure 3.2). 
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4.0 DISCUSSION OF OCEAN BACKSCATTERED SIGNAL CHARACTERISTICS 

This section presents a discussion and interpretation of GEOS-3 data for 
cases in which the AGC values increase markedly for brief periods. These 
abrupt AGC changes have been observed in a considerable number of passes; 
Figure 4.1 shows the geographical distribution of occurrence of these AGC 
changes noted in the examination of approximately 75 records. It was thought 
that such a map might show a pattern, or grouping, in these occurrences (such 
as near the Gulf Stream), however, the distribution shown in Figure 4.1 is 
considered to essentially represent the geographic distribution of the under- 
lying data base. This result suggests that these comparatively brief eleva- ~ 
tions in AGC values are the result of relatively calm ocean surface conditions, 
and are not caused by anomalous conditions. This premise is examined in the 
following paragraphs. 

First examining the theoretically predicted values of a?; standard ref- 
erences give the appropriate form as 

- ,„«->= 

s 

where R ® Fresnel reflection coefficient, 

l - surface height correlation length, 

a ~ rms surface height, and 
s 

0 = off-nadir angle. 

2 

For the GEOS-3 case: -cos0 - 1, sin 0 - cr/h and assuming R - 1 



r- 2 

x tan 0 

40 
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■ Tabulated values for this equation are given below for x = 60 ns. 


z_ 

a 

s 

a (dB) 
o 

3 

3.52 

5 

7.96 

10 

13.97 

50 

27.9 

100 

33.75 

500 

42.2 


The above value of t was used to provide numerical values applicable to 
the AGC gate region. Note that the angular displacement of the AGC gate lo- 

cation from the nadir point is ~ 0.27 degrees, for which the decay effect of 

✓ 

the antenna pattern is ~0.1 dB. The above tabulated b° values indicate the 
values of a 0 observed during the AGC step-changes are entirely consistent 
with backscatter theory. These theoretical values indicate that the elevated 

■# 

AGC conditions correspond to very calm surface conditions or to swell-dominated 
seas. 

Figure 4.2 shows computed waveshapes as a function of Z/<y . These wave- ' 

s 

shapes show the expected decay in the plateau region between the leading edge 

and the AGC gate, as a function of the ratio of surface correlation length (Z) 

divided by rms surface roughness (a ) . This result suggests the possibility 

s 

of using attitude/specular gate and data as a means of estimating domi- 

nant surface wavelength. This possibility is next investigated. 

4.1 Attitude/ Specular Gate Behavior 

The change in the -attitude/specular gate value as a function of Z/o 

s 

may be computed using the previously given equation for 0° and the known an- 
tenna pattern behavior. Taking 0 in this equation as ~1.0 degree for the 
attitude/specular gate angular location, the antenna pattern effect to be 3dB 
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RELATIVE VOLTAGE 



■ Figure 4.2. Computed GEOS-3 video waveshapes as .a function of the ratio of 
ocean-surface' correlation Length to rms wavehei'ght. ; 


and the AGO gate nominal value to be 100 rav, the following attitude/specular 
(A/S) gate values are obtained. 


l/o 

s 

A/ S (rav) 

20 

48 

40 

45 

80 

31 

100 

23 

200 

2 


In order to calculate correlation length as a function of attitude/ 
specular gate voltage, V , a polynomial was curve-fitted to the above 

3.S 

values. The results were 

H = H. [*54.33 -2.27 V +.0587 V 2 ' - 6.79 * 10“* V 3 1 
1/3 L as as. as J 

where H, 4(X is in meters and V .is in rav. 

1/3 s as 

When this procedure was implemented and wavelength computed using mea- 
sured V values it was soon found that GEOS-3 attitude/specular (A/S) data 

flS 

shows the following pattern: V values tend to first increase above 50 mv 

as 

when an AGC step-change is encountered, to decrease in an expected fashion, 
and to again increase above 50 mv as the AGC feature is exited. This char- 
acteristic is shown in the data given in Table I, which was recorded over the 
Gulf of Mexico, orbit 1164. 

Since backscatter theory does not admit to an increase in energy with 
off-nadir angle, this behavior was thought to be caused by the geometry of 
the A/S gate. Figure 4.3 shows the relative spatial areas of the sampling 
gates. Note that along-track signal level changes can be anticipated by the 
A/S gate by ~2 seconds (or -1 low data rate frame period). A computer pro- 
gram was written to simulate the effect of such AGC changes on the A/S gate 
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TABLE I 


Ln 

N> 



Q. wj 
& 'y 


LATITUDE 

Degrees 

LONGITUDE 
’ Degrees 

a ■ 

’ n , 

Meters 

AGC 

dBm 

H l/3 

Meters 

AASG 

m.v. 

26.162 

277.793 

• 47 

-70.80 

-2.95 

49.2 

‘26.055 

277.722 

.53 

-70.69 

-1.78 

48.4 

25.948 

277.650 

.60 

-70.59 

1.84 

50.3 

25.841 

277.579 

.52 

-70.15 

-2.52 

52.8 

25 . 735 

277.508 

.35 

-69.09 

.47 

48.8 

25.628 

277.437 

.47 

-69.28 

-2.67 

50.3 

25.521 

277.366 

.99 

-69.21 

-1.39 

54.7 

25.414 

277.295 

.83 

-67.89 

2.14 

48.8 , 

25.307 

277.224 

.43 

-68.42 

1.38 

50.3 

■25.200 

277.153 

.42 

-68.22 

-3.34 

54.6 

25.092 

277.083 

.41 

-68.26 

-3.63 

58.3 

24.985 

277.012 

.52 

-66.42 

-2.90 

53.9 

24.878 

276.942 

.37 

-66.97 

-3.85 

49.5 

24.771 

276.872 

.82 

-67.63 

i. 45 

62.0 

24. 664 

276.802 

.55 

-67.30 

-1.75 

75.1 

24.557 

276.732 

.54 

-66.05 

.3.27 

79.6 

24.449 

276.662 

1.06 

-63.09 

3.40 

77.0 

24.342 

276.592 

.64 

-60.76 

3.04 

71.1 

24.235 

276.522 

.68 

-57.84 

-1. 52 

44.8 

24.128 

276.452 

.4‘4 

-58.26 

^3.46 

44.8 

24.020 

276.383 

.60 

-57.41 

-2.70 

40.8 

23.913 

276.314 

.38 

-58.17 

-3.28 : 

36. 4‘ 

23.805 

276.244 

.61 

-59.20 

-3.88 

43.0 

'23.698 

276.175 

.56 

-58.68 

-3.61 

46.3 

23.590 

276.106 

.85 

-59:53 

-3.45 

40.1 

23 . 483 

276.037 

.57 

-62', 85 

-3.03 

63.8 

.23.376 

275.968 

.47 

-61.69 

-1.54 

55.0 

23.268 

275.899 

.71 

-62.65 

3.94 

49.6 

23.160 1 

275.831 

.62 

-64.35 

-1.44 

45.5 

23.053 ■ 

275.762 

.62 

-65.60 

-2.28 

47.0 

22.945 

275.693 

.88 

-67.54 

-2.03 

47.8 

22 . 838 

275.625 

.39 

-68.33 

-2.54 

50.3 

22 . 730 ' 

275.557 

.58 

-69.24 

-3.43 

50. 7 

22.622 - 

275.488 

.58 

-69.64 

-1.92 

51.7 

22.515 - 

275.420 

.6-2 

-69.58 

-1.87 

51.0 

22.407 

275.352 

.87 

-70.12 

1.70 

52.1 

22.299 

275.284 

.28 

-69.79 

-1.35 

51.0 
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Figure 4.3. R^ - inner radius and = outer radius of footprint. 
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behavior, for the case of ocean surface conditions which are homogeneous in 
the across track direction. This program comprised the following operations: 


1) AGC values a, , a„ , . . . a were read in 

1 2 n 

2) these dBm values were converted to power values A 

, n 

3) the simulated A/S level was calculated based on the AGC 
power values one frame ahead and behind the nadir point, i.e. 


_ A(n-l) + A(n+1) 


4) these calculated A/S values were normalized to the nadir 

t 

value and to 50 m.v. using 


A(n) 

A/S 


x 50*' 


\ 


A typical result from this simple model’ is shown in Figure 4.4. The low- ; 

> - ' * ' t 

er graph in this figure shows the AGC features which were, recorded in the ' 

£ 

Mediterranean Sea near the island of Crete (orbit 3469 ),r Note that there are . 

- t 

two brief changes of ~ 10 dB in AGC value which" are not of sufficient dura- / 

J 

/ 

tion to effect non- transitory conditions in the AGC and A/S gate footprints.' 

j 

i t 

The upper graphs in this figure show both calculated and measured A/S gats’ 

responses. These curves show that AGC disturbances alone can cause eleva- 

* * s 

tions above 50 m.v. in the A/S gate value. Although the calculated and observ- 
ed A/S gate values are not in close agreement, the peaks and troughs and over- 
all characteristics are considered to be in sufficient correspondence to vali- 
date the postulated mechanism. The. lack of precise agreement, is attributed 
to the (unknown) across-track variations in ocean ba.ckscatter . 

I 

The above results demonstrate that the analysis of waveform variations 
inferred by A/S gate-data requires that the ocean surface be essentially 
homogeneous over a spatial extent of tens-of-kilometers . Figure 4.5 shows 
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AGC (dBm) 



Figure 4.4. Comparison of AGC values with measured and calculated 
Attitude/Specular gate values, rev. 3469. 
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' Figure 4.5. An example of a protracted rise in AGC level'.' 


A AS 6 VOLTAGE ( M.V.) 


cr° and A/S gate data for one such record; these data were recorded over water 
between Cuba and Florida. Note that the elevated <J° values in this record 
persist for ~ 30 seconds. Measured significant waveheight during this period 
was *- 1 meter; use of the above equation fox correlation length gives £ = 33.7 
meters. Note that this value coupled with the high Cf° value argues that 

swell is dominant; also the observed <7° value is in very good agreement with 
the theoretical value for this £/cr ratio. Using the approximation: wave- 

length X = 4& and the standard equation for wave period T, 


T = 



yields X = 9 . 3 seconds , which is a reasonable period for swell conditions . 

The main uncertainty in the above calculations is in the ^ value. A 
total of 10 per-frame H^y ^ values were averaged to obtain the value quoted 
above. The standard deviation in these 10 per-frame H^y^ va ^- ues was 
meters. Based on these results, it is concluded that the present H^y ^ al- 
gorithm is probably not sufficiently stable (at low waveheights) for use in 
estimating surface correlation length. The variance-based H^y ^ algorithm, 
discussed elsewhere in this report, may prove to be adequate for this purpose, 
if surface correlation length estimates prove to be a useful parameter. 


4.2 Spatial Variability of H l/3 Estimates 

Figure 4.6 illustrates the autocorrelation function of frame averaged 
for two different averaging times. When the frame averaging time is 

, * t 

changed from five to three the decorrelation time is affected by more than 
60 percent. This is a significant effect and indicates that a five frame 
average of values can significantly affect the statistics of a three 

frame averaged R. process. 
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Figure 4.7 illustrates the effects of a five-frame average of on 

the autocorrelation function. The solid curve for the autocorrelation func- 
tion of unaveraged data clearly indicates the presence in the data of 

two uncorrelated processes. The first of these processes rapidly decorre- 
lates and is probably due to Rayleigh noise caused by return signal fluctua- 
tion. The second process is characterized by slow decorrelation and repre- 
sents sea state effects. Note, however, that in the average curve the two 
processes have been smeared together so that they are essentially indistin- 
guishable. Based on these results it is concluded that care is required in 
interpretation of results based upon filtered data. The results of this 

paragraph indicate that can be significantly altered by filtering opera- 

tions . 

4. 3 AGC Calibration 

This paragraph gives results of values computed using both the -'clean 
signal" and the "clutter signal" AGC calibrations. Figure 4.8 shows C a values 
obtained under moderate waveheight conditions and Figure 4.9 shows 0° values 
during a period in which the AGC values experience a brief step change of 
-10 dB. As shown in Figure 4.1, such changes have been observed over a rather 
wide range of geographic locations. Figures 4.8 and 4.9 show the clutter 
calibration data to yield better agreement with published 0° values than do 
the clean calibrations. For this reason, it is recommended that "clutter 

v 

signal" calibrations be used in all a° computations. 
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Figure 4.8. Radar cross-section values for moderate sea states. (~6 meters, H.^) 
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5.0 ALTIMETER DATA PROCESSING CONSIDERATIONS 


Experience with the GEOS-3 data has shown that preprocessing or data 
editing is required before optimal filtering procedures can be applied. Auto- 
mated techniques have been developed which provide for wild-point interpola- 
tion and correction and for estimation of residual orbit bias and slope pa- 
rameters [1]. These parameters are obtained using statistical regression- 
techniques. 

Non linear characteristics of the GEOS-3 split-gate tracker have been 
found to cause negative asymmetries in the altitude data. Figure 5.1 shows 
altitude data for orbit 2023 over the- Gulf of Mexico; this record displays 
the characteristic negative (or downward) perturbations in the altitude data. 
Figure 5.2 reproduces a section of the orbit 2023 data on an expanded scale 
and containing threshold data. The threshold altitude algorithm is seen to 
essentially remove the altitude disturbances and to indicate that they are 
hardware-specific in origin. Based on .this behavior, and because low data 
rate telemetry data cannot use the threshold algorithm, computerized edit 
procedures were implemented to interpolate through such periods. 

5.1 Sea State Bias 

The term sea state bias has been used to account for any systematic dif- 
ferences between mean value of the geometrical ocean surface and radar sensed 
mean value. This type of bias is thought to arise due to either the trochoi- 
dal shape of ocean waves or the increased occurrence of capillaries on wave 
peaks compared to wave troughs. Either of these factors could cause an in- 
crease in radar cross section per unit area with increasing distance below 

* 

the wave crests. The direct measurement of this bias term under deep water, 
long fetch conditions would be extremely difficult. It would require the ac- 
quisition of simultaneous waveheight and radar backscatter data of very high 
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absolute and temporal accuracy in mid ocean, under high sea state conditions. 
The only known instance in which direct measurement has been attempted is the 
work of Yaplee, et ?al. , reported in 1971 [2]. These measurements, conducted 
from a tower located off the coast near Norfolk, Va. , showed a bias toward 
the wave troughs of magnitude 

At - .17 H 1/3 - .31 

where At is the time bias in nanoseconds and is the significant wave- 

height in feet. The empirical relationship shows a bias change of ~ 8 cm 
for a waveheight change of ’1 meter. This degree of bias will be shown to be 
in general agreement with the bias results to be discussed. 

Tn contrast with the results reported in [2], the sea state bias results 
to be given below relate to instrument-induced biases as well as ocean sur- 
face effects. These instrument biases are primarily caused by changes in the 
tracker equilibrium point which are in-turn induced by any changes in the 
received waveshape. Such waveshape changes are caused by variations in: 

■ a. received signal-to-noise ratio, and 

- b. significant waveheight. 

*■ 

- ' ' i ^ 

Experimentally measured values of altitude bias as a function of signal-to- 

noise ratio are in the range of 2.0 ns, see [3]. Bias effects due to pointing 

> 

angle and significant waveheight cause corresponding changes in the split 
-gate tracker equilibrium point; these effects have been analyzed in [4]. 

The sea state biases reported here are based on analysis of data from 
near-overlapping satellite passes. Altimetric data for such passes is first 
examined and, if necessary, edited to remove data discontinuities and then 
filtered using a Wiener convolutional algorithm [5]. The ensuing data from 
several passes is then adjusted in absolute altitude .value to provide over- 
lap for the low sea-state portions of the passes.,- Figure 5.-3 shows raw 
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Figure 5.3* Unprocessed sea surface height data for two near-over lap ping passes* 



altitude data (10/sec) for two passes; orbit numbers 1633 and 6893 for a 
» 

North Atlantic segment over the latitude and longitude limits shown on the 
abscissa. Orbit 6893 data comprised the high sea state pass and the data is 
considerably noisier than orbit 1633 data; however, • the permanent geoidal 
characteristics are seen to exhibit very good repeatability. Figure 5.4 
shows both the filtered sea surface'; height and the sea' state data relevant to 
these data periods. Note in the lower figure that the'.'passes (which have been 
offset by - 2 'meters) show a systematic departure which correlates with the 
high sea state period shown in the upper figure. „ -3 

The above analysis technique has been used with -a total of five orbits 
to produce the results shown in Figure 5.5. These results show the altitude 
measurement to be increased by the presence of high seas; that is, the measur— 

v 

ed mean- sea-level value is depressed downward due. to high seas. Based on the 
results shown in Figure 5.5, the bias effect increases with increasing sig- 
nificant waveheight and i s in the order of 10 percent of the 

value for 10 meter seas. 

5.2 Precipitation Sensitivity 

Figure 5.6 gives calculated values of 1) single pulse signal levels and 
2) signal- to-noise ratio for one second post detection averaging, both as a 
function of precipitation level. These results indicate that moderate to 
—heavy- precipitation shouid'be" detectableTin - the Global mode. Other consid- 
erations, however, indicate that hardware changes would be necessary to im- 
plement this capability; the present noise gate is an a-c coupled circuit 
that does not provide a direct measure of noise power. In the event that the 
GEOS-3 backup system is used,, simple changes in the noise or sampling gates 
and the AGC loop would permit observation of precipitation return. 
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Figure 5.6 » Calculated G~Mode and 1— Mode received SNR versus rainfall 
rate for one second post detection integration „ 



5.3 Waveform Sampler Corrections 

The following paragraphs summarize our efforts toward deriving correc- 

tions for the sampling gates. The initial approach to this problem’ was to- 

* \ * - „ f 

use BIT/CAL data to identify offset and gain differences between samplers. 

As discussed below, this approach was not feasible because BIT/CAL offset 

values did not correlate with data acquisition offset values. The bias-test 

, *( 

* Y * „ 

portion of the BIT/CAL data also showed waveform changes, especially in the 
12 — gate L values. It is not known whether these were due to sampler effects 
or changes in transmitter leakage. 

v 

The second approach comprised the extraction of offset corrections from 
fitting an error function to the waveform and gain corrections from the stan- 
dard-deviation waveform. Two difficulties were experienced in this approach: 
a- study of gate-to-gate covariance indicated that' timing differences were 
also present.' In addition, offset corrections derived from the on-board 
averaged waveforms and the software averaged waveforms were not in agreement. 
At this point it was decided that gain corrections could not ’be obtained for 
the low* data rate case; and the approach taken was to obtain offset correc- 
tions from the error function fit and timing corrections from the covariance 
calculations. Shortly after this approach was taken it was learned that E. J. 
Walsh of WFC had obtained pre-launch timing corrections from G. E. and had 

derived offset corrections, both of which generally agreed with our results. 

< . 

Our activities were discontinued at this point . 

As the results given below show, empirically derived offset corrections 
appear to be valid to within -2 m.v. The 2 m.v. uncertainty derives from 
short term drifts in the offset values. ' . 

Figure 5.7 shows graphed values of BIT/CAL -and backscafctered waveforms 
for Pass 217. Data from other passes have also been examined and the BIT/CAL 
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Figure 5.7. Comparison of BXT/CAL and experimental 
waveforms; Pass 217 (6-Frame average). 
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waveforms were found to be similar to those shown in Figure 5.7. In this 
figure the baseline is displaced upward one division for each of the wave- 
shapes shown. The following comments pertain to the data shown in Figure 5.7 

1) The unevenness of the ocean scattered waveshape is much 
larger than the expected statistical uncertainty. The 
6-frame averages contain approximately 1900 samples , there- 

- fore, the one-sigma confidence bound should be approximately 

2 %. 

2) There .is no clear cut relationship between the noise region 
sample values for the Bias Test and the experimental waveform. 

3) The fact that the IF and video waveforms are locally uneven 
and the differences between the voltage and standard devia- 
tion relative waveshapes suggest that dc biases and gain 
d i fferences both exist . 

Figure 5.8 shows a time-sequence graph of several "noise" or early gates 
these results partially explain the lack of agreement in bias values noted 
above. The upper curve shows per-frame values for ARS-4 , which represents 
the- on-board averaged values for the fourth gate-. The next three curves 
show the software averaged values (AW) for gates 1, 2 and 4. These three 
curves are rather - surprising for two reasons: The apparent correlation be- 

tween gates 2 and 4 is striking; and the correlation- with- gate 1 i's higher' 
than expected. The rate of fluctuation is also much greater than expected; 
if these off-sets were due to d-c drift, temperature, etc., the fluctuations 
would be much slower. The lower graph shows AGC values and the largest gate 
excursion is seen to coincide with the highest AGC value; thus the noise gate 
values are seen to increase during a period in which they would be expected 
to show a decrease. The level of correlation present in the gate values 
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FRAME NUMBER 
(ORBIT 2762) 


Figure 5.8. Comparison of on-board ARS and software AW averaged gate values with AGC data. 

(For 320 values of intra-frame data, tape dump showed Gate 2 rms of 4 m.v.) 



indicate that the changes are due to system variations rather than signal 
information. The per-pulse values for one frame were examined and quantiza- 

rw 1 1 ' ' 

• i < 

tion effects were strongly in evidence. The data showed typically 3 step or 
transition values of ~ 7 m.v. * Inspection of the variances of the per-frame 
values shows that -the quantization distribution is not uniform. Thus the 
step changes appear to be a slowly-varying effect. 

t 

Figure 5.9 shows the computed covariance for the first five gates based 

on use of the per-pulse values. The signal variability in these gates are 

* § ; 

largely due to offset drift as; depicted by the ~ 7 m.v. quantization. This 

1 

figure shows gates 2 and 3 to -be much more highly correlated than gates 1 and 

2. This implies a wider gate separation in the latter case. Figure 5.10 

* 

shows the covariance in the plateau region for- two cases. Since gate 14 was 
selected as the reference gate,'. the right-left asymmetry in' the covariance 

s 

function is indicative of sample timing differences.. 

A set of gate ^spacings were derived by least-squares fitting a Gaussian 
. _/ 
covariance model to the experimental data; the fit was constrained to re- 
quire that the sum of the gate spacings approximate the design value (4x6; 25 
n.s.J. The fit obtained is shown in Figure 5 .11 and the derived gate spacings 
were ’ ^ 

Gate No. 12 13 14 !L5 ‘ 16 

_Spacing_ (n. s^.). _ _ _ _5. - — 9- - — 5- — - 6 -. 5 - ------- - 

* 1 

With the exception of the gates 12-13 spacing, these results are in good 
agreement with Walsh's results (Figure 5ll2). 

5.4 Tracking Jitter Correlation Properties 

Autocorrelation and tracking jitter probability density functions were 
extracted from the orbit.. 578 data. ■■ Figure 5.13 shows the pre-launch measured 
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Figure 5.11. Least-squares fit of Gaussian 
curve to experimental data. 
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Figure 5.12. Comparison of offset corrected waveforms with, uncorrected waveforms. (From Walsh) 
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tracking jitter correlation coefficient. Figure 5.14 shows (1) the computed 
correlation for orbit 578 and (2) the over land correlation for orbit 2236. 

The agreement between --the orbit 578 data and the pre-launch data is quite 

? 

good; this over-ocean- experimental result shows the loop bandwidth to be in 
the neighborhood of 5 Hz, based on the correlation period shown. The over 
land result for this data-span shows that the very different mean waveshapes 
encountered cause the time discriminator characteristic and thus the loop 
bandwidth to differ substantially from the over ocean values. It should be 
noted that waveform data for the overland segment of 2236 [2] shows consider- 
able saturation and the tracking loop may be operating non-linearly. 
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Figure 5.14. Computed correlation coefficients of altitude residuals. 


6.0 CONCLUSIONS AND RECOMMENDATIONS 

Comparisons between the deconvolution sea state algorithm and the WFC 
algorithm show good agreement and it is concluded that the WFC algorithm is 
entirely adequate for GEOS-3 sea state computation. Use of the deconvolu- 

i , ! 

tion algorithm has,. shown that extraction of the raclar inferred, waveheight 

* * ^ > 

density, function from GEOS-3 data is possible under 'high sea state conditions. 
These densities suggest that shorter pulse, higher. data-rate sensors (such 
as the SEA SAT altimeter) will be able to resolve the' underlying. densities 
in the low to moderate sea, state range. It is recommended that the. .decon- 
volution algorithm given herein be utilized with SEA SAT data to investigate 
the possible improvements in accuracy of sea state estimation compared to al- 
gorithms that assume that the received waveshape is strictly an error func- 
tion form. (Time sidelobe effects on the SEA SAT data should also be con- 
sidered.) 

Both the Cramer Rao error analysis and examination of sea state spatial 
variability have shown that the achievable sea state resolution is bounded 
by either Rayleigh noise or by ocean surface inhomogeneity. An averaging 
period of , ~ 20 seconds is considered to be near— optimum for GEOS-3 data pro- 
cessing. 

Comparisons between theoretical a° values and GEOS-3 observed values 
have been found to-be highly consistent. The "step changes" sometimes ob- 
served in o° data, in- areas such as the Gulf Stream, were found to corre— 

* 2 

spond to theoretical values for which the.,rms surface slope is small. Since 
surface slope is dominated by waves in the capillary range, such calculations 
indicate that these cr° changes will occur over either very calm or swell domi- 
nated seas. Based on! the comparisons given, it is recommended that the "clut- 

i *,*■•••* *'■*••* * " * 

ter" AGC calibration -be used in calculating cr° values. 
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Data analyses reported herein show a systematic change, . or bias, in in- 
dicated sea surface height as sea state increases. This change is generally 
less than 10 percent of the significant waveheight and, as such, is consider- 
ably less than orbit errors. Should orbit errors dramatically improve in 
future systems, this error source should be re-examined. 

The presence of d-c biases in the waveform samplers have been found to be 
a source of difficulty in waveform data processing. Especially frustrating 
is the apparent inconsistencies between the calibrate and data acquisition 
indicated biases. In future systems, it is recommended that laboratory studies 
be conducted to resolve these differences arid that use of chopper stabilized 
circuits be considered. 

The GE0S-3 system has proven to be a reliable and versatile satellite 
system and sensor; consideration should be given to possible' use. of - the back- 
up hardware system in other remote sensing investigations. Some of the poten- 
tial uses are as follows: 

1. With modifications to the pulse length and repetition rate ( ~6 ns 
and 500 per sec.) the sea state resolution. could be increased by a 
factor of 5 . With these comparatively inexpensive changes the sys- 
tem could be orbited in conjunction with other Shuttle launches to 
provide a denser grid of sea state and surface wind data compared 

to that available from the SEA. SAT altimeter. Experience with GEOS-3 
has shown that the principal limitation to its sea state measurement, 
capability is in the coverage provided by a single (nadir) sensor. 

2. With primarily altitude tracker and AGC modifications the system 
could be optimized for data acquisition over ice and land areas . 

Other studies [1] have discussed the usefulness of global ice boundary 
data and ice sheet topography measurement to NASA climate programs. 
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3. The addition of a multiple feed structure and possibly a larger an- 
tenna would provide wind data over perhaps ±100 km off nadir. This ' 
would yield useful information on surface wind- homogeneity and, to 
some degree, wind direction. 

4. All of the above sensor concepts pertain to capabilities-- of the 

Intensive Mode. Analyses summarized in this report indicate that 
minor 1 modification of the Global Mode would permit the acquisition 
of' precipitation data. This capability would thus be incorporated-' ■ 
with any- of the above system modifications. ■■ - 

In addition to these hardware-specific suggestions, the GEOS-3 program, in--'-' 
proving the topographic and sea state measurement capabilities- of- short pulse"- 
al time try, ‘ has- demonstrated the greatest deficiency to be that of ■ area ; cover- 
age. Since finite-swath altimeter concepts have' not been forthcoming (for--- 

✓ 

practical antenna sizes) , the obvious solution appears to be use of multiple- ■ 
satellites. For sea state measurement only, consideration - should be given- to'- 
sensing techniques other than those requiring • short pulse-lengths. Additional 
studies are needed in these areas. 

REFERENCES 

[1] Miller, L;-S., Brown, G. S. and Priester* R. W. : "A Summary of GEOS-3-"' . 

— Data- -Studies^ Applied' 'Science Associates , Inc., Technical Report, Novem- 
ber, 1977. 
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Appendix A 


Computer Program Listing 

The results presented in paragraph 2.3 were obtained using the Fortran 
program which will be described in this appendix. In the attached listing 
Fortran statements often appear to begin in position 9 rather than the stan- 
dard position 7 . This is a result of the TAB key feature supported by the 
R.T11 TEXT EDITOR in DEC’s LSI 11/03 computer. If storage space is. .a problem, 
an overlay structure can be used as illustrated- in Figure A.l. Such a struc- 
ture has been used under RT-11 for which TTY I/O is designated by devices 
5 (input) fl (output) respectively in the listings contained herein. The list- 
ing presented in this appendix was used to obtain the results of paragraph 
2.4.1; therefore, the attached listing can be used in its present form to 
obtain' anc * Mt) if H^y^ 4 m. To use the code as described in para- 
graph 1.4.2 two changes are required: 

(a) In subroutine H13MSR the instruction DATA Ul/30*1.0 

should be used in place of DATA Ul/0, .25, .5, .75, 26*1./ 

(b) In- subroutine WAVEHT change statement marked with * 
(immediately after statement label 90) to read as follows: 

SWH=DABS (SWH) *SWH -57.0 

An example of executing the code immediately follows the FORTRAN listing 
and will be described at this time. The input command R H13TST causes the 
executable code corresponding to the source program to be loaded. The pro- 
gram then requests that NNUM , NDENOM , NINP be input in format 313. NNUM = 
m-1 (where m = order of H(s) numerator polynomial), NDENOM = order of H(s) 
denominator polynomial), and NINP is an integer which specifies the number of 
sample gate values which will be input (usually this will be 16 for GE0S-3) . 
Normalized gate samples are printed next. These values are computed from the 
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Overlay structure of Fortran subroutines 
used to obtain results of paragraph 1.4’i 



gate samples input earlier as follows: 

(a) The average of the first four samples is subtracted 
from all gate values in order to remove offset effects. 

(b) The input set of samples is expanded in number to 25 
by padding from NINP to 25. The value inserted is 
specified by the average of input samples NINP , NINP-1, 
and NINP-2. This is done in order to enhance the sta- 
bility of the H(s) identified by the algorithm. 

After the identification algorithm has executed successfully, the pa- 
rameter vector is printed. The parameters are listed in the following order: 
a l* a 2’ **’ a n * ~^l s “^ 2 * * * * ,_ ^n‘ ^ ext probability density is output; the 
numbers printed must be divided by 10. Significant waveheight (SWH) , m2(AT2), 
m^(AT), and A (a) are printed on the next line (see paragraph 2.3 for defini- 
tion of these quantities). 

Residuals are defined as input gate values minus model output resulting 
from identification procedure. The final output is the sum squared error, 
SSQER. 
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C 

C MAIN PROGRAM 

C 

DIMENSION ZC64) jTR<64>*T HC64) 

DIMENSION COVC 144),C0VI C 144) 

DIMENSION XH<350),XHT<350> 

DOUBLE PRECISION 7l> TR> TH, CO V> COVI > XH, XHT 
WKITEC 7*501) 'I 

501 FORMAT* ’ INPUT NNUM* NDENOM»NINP-FOHMAT 313') 

READ* 5* 500) NNUMER* NDENOM * N I NP 
N DATA- 2 5 

500 FORMAT (31 3) 

C NNUNER=ORDER OF N(S)=ZER0 IF ONLY CONSTANT IN NUMERATOR'. 

C NDEN0M=0 RDER OF DENOMINATOR* MAY BE 2>3*4*0H 5 ONLY 

C NDATa=NO. OF DATA POINTS TO BE PROCESSED 

C NINP=N0 » OF G EOS-3 INPUT SAMPLES 

NN-N.N UMER 
ND=N,DENOM 
NDAT=NDATA 
NPR=,NNUMER+N DENOM+ 1 
NPR=NPR+0 

C Z*TR*TH*XH*XHT*COV*COVI. ARE STORAGE AREAS F0R'H13MSR 

CALL. H 1 3MSRC Z* TR* TH* XH> XHT* COV* COVI * NN*ND*NDAT*NPR*NINP) 

STOP, 

END i ‘ • 

SUBROUTINE H 1 3MSh( ZETA>THMAT > THETA* X-HAT * XHATT > COU* CO VI > NNUMER* NDEN r - 
COM* N DATA* npram* NINP > 

C ZETA=MATRIX: (A** (-!)*( THETA-TSTEP*I ) ) 

C THETA=M ATRIX : (A**(-l))*( PHI - 1 ) 

C PHI = t'FiM AT = STATE TRANSITION MATRIX 

C XHAT=MATRI X OF PREDICTED STATE VECTORS 

C XHATT=MATRIX TRANSPOSE OF XHAT 

C THIS' PROGRAM ESTIMATES Hl/3 FROM 16 GATE SAMPLES 
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DIMENSION PI (8>4)*YOC30).» YM (30>.»U(30> 

DIMENSION X1K(83,X1KP1 (8 )> X2K( 8 )> X2KP1 <8> 

DIMENSION PBACKC8) »XI I < 30 ) > BI GY C 30 > , ERV< 30 > , P ( 1 2 > , YSOL< 1 2 > 
•DIMENSION COVCNPRAMjNPRAM ) » CO VI (NPRAM..NPRAM ) 

DIMENSION THETACNDENQM* NDENOM ) > ZETA<NDENOM# NDENOM ) 

DIMENSION TRMAT (NDENOM..NDENOM) 

DIMENSION D1 C64),,D2<64>, D3C64) 

DIMENSION Y C 30 )>U1<30) 

DIMENSION PROBC30) 

D1,D2,D3 ARE TEMP. STORAGE FOR MATRICES COMPUTED IN SUB-ROUTINES 
THETAS USES C=D1,AINV=D2 
MATEXP USES DMeDl , FM=D2, XMAT= D3 

DIMENSION XH AT C N D AT A., N P R AM > > X H AT T < N P R AM N D AT A ) 

DOUBLE PRECISION ZETA, TRMAT/tHETA, XHAT, XHATT, COV, COVI 
DOUBLE PRECISION YOj YMj U^X 1 K.»X 1KP 1 > X2K..X2KP I » PBACK> XI I 
DOUBLE PRECISION BI GY* ERV, P, YSOL, D1 , D2> D3* YBAR> SSQER 
DOUBLE PRECISION SSQERBj TSTEP> EPS 1 > EPS 3* FAG 
DOUBLE PRECISION ERCHK* DU* DY 
DOUBLE PRECISION XEXP>YHI 
DOUBLE PRECISION SWLAT2,AT,A 
SET UP INPUT FUNCTION 

DATA Ul/0.* .0288, . 195.5 j . 5000* . 8045, . 97 1 2, 24* 1 . / 

DATA Ul/30*1.0/ 

DATA Ul/0« j »25.» *5j*75j26*1»/ 

C ' DEFINE INITIAL ESTIMATES OF PARAMETER VECTOR 

DAT A Pl/3»>3.^6*0.^ 15 . » 15* > 6 . * 5 * 0 • * 1 0 5 • .» 105.»45.* 10. ..4*0* >945.^945 
+ »j420«j 105* » 1 5 • .» 3*0 . / • 

C READ GEOS GATE SAMPLES 

WRI TE ( 7 f 500) 

500 FORMATC-* INPUT GATE SAMPLES ’/5< * . ’)) 

- READ(5^501)(Y(I)>I=UNINP) 

501 FORMAT C 5F 1 0.5) 

DO 5 I = 1 jNDATA 
U ( I > = DBL E C U 1 < I ) ) 

DO 6 I=1,NINP 


5 
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i. 


6 

C 


C 


10 


1 1 
S00 
801 


C 


12 

C 


Y0( I ) = DBLECYCI-> > 

PREPRGCES'S INPUT DATA SAMPLES 
YBAR=VO< i >+YOC2)+YO<3>+YO<4> 

YBAR=YBAR/4. 

SUBTRACT OUT- YEAR 
DO 10 1= 1 jNINP 
YOU ) =YO < I ) -YBAR 
CONTINUE 

YHl =YO (NI MP ) +YO CN I NP- 1 > +YO (N INPtS > 

YHI=YHI/3» 

NINP1=NINP+1 ' ‘ 

DO 11 I=NINF1,NDATA 

YOCI)=YHI 

WRI TEC 7 #800) 

FOftMAT'C ’ NORMALIZED GATE SAMPLES 1 
WRI TEC 7 .>801) ( YO ( I > , I = 1 * NDATA) 

FORMAT < 8F8 • A ) 

YO CNDATA+1 )=YO (NDATA) 

YO VALUES NOW REPRESENT UN-NORMALIZED PROB - DISTRIBUTION 
SSQERB=0. 

TSTEP=', 1 
NMODE=£0 • 

EPS 1= 1 * D- 1 3 
EPS 3= 1 • D- 1 3 

APPLY SCALINGS TO SELECTED P-VECTOR AND INITIALIZE ESTIMATE 
VECTOR' 

FAC=2.\ 

DO 12 1 = 1 /NDENOM'- 
I EX‘P=NDENOM+ 1 - I . 

INDX=NNUMEK+1+I 
ININDX=NDENOM- 1 

PC INDX):=DBLE(P 1 C I j ININDX) ) *( (FAC ) **I EXP> 

P C I N DX )’ = - P ( I N DX ) - 
IFCI.LE.NNUMER+1) P(I') = 0:. ■ 

CONTINUE • ’ ' 

INITIALIZE BACK-VALUED P-VECTOR 
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-J 


36 

C 


15 

C 


C 


C 


25 

30 


C 

c 

c 

c 

c 

C 35 
C 


C 


C 


DO 36 I=l,NDENOM 
1 1 = I +NNUMER+ 1 
PBACKCI >=PCII ) 

CONTINUE 

TEST FOR MONOTICITY OF YO 


NOSW=0 


DO 15 I=6,NDATA 

I F( YO ( I > *GT • YO ( I + 1 ) ) NOSW=NOSW+l 
CONTINUE 


NO SW=NO • OF TIMES THAT CURVE WAS DECREASING 


ITEHNO=0 

********* TH I S IS MAIN LOOP POINT********** 

DO 150 ILOOP=-1,16 
I TERNO=I TERNO+1 

DECIDE IF MODEL CYM) OUTPUT OR EXPERIMENT (YO) OUTPUT SHOULD 

BE USED IN COMPUTATIONS 

IF(ITERNO.GE.NMODE) GO TO 30 

DON'T USE MODEL OUTPUT YET 

N DAT P 1 =N DATA* 1 

DO 25 1 = 1 , NDATP 1 

YMCI )=YO (I ) 

CONTINUE 

CONTINUE 

DO 35 1 = 1 , NDENOM 
X 1K( I )=0» 

X1KP1 Cl )=0. 

X2KCI >=0. 

X2KP1 (I >=0. 

CONTINUE 

COMPUTE TRANSITION MATRIX 

CALL MATEXPCTRMAT, P, TSTEP,NDENOM,NNUMER,NPRAM, D 1, D2, D3> 
COMPUTE THETA FOR HCS> SUB- I 

CALL THETAS ( THETA, P, TRMAT, NDENOM, NNUMERi NPRANJ, D 1 , D2 > 

COMPUTE ZETA FOR HCS) SUB-I 

CALL Z ETAS (ZETA, P, THETA, TSTEP, NDENOM, NNUMER,NPRAM, D1,D2> 

JXK= 1 


o c 

8 S' 

JO 

cj r? 

o 

tej 

8 


P 



r. 

oo 


C SKIP • PROPAGATION AFTER 1ST ’ITERATION 

IFCITERNO.GT. 1 ) GO TO 39 

C NOW C.ALL VECPRP TO PROPAGATE STATE AND FORM XHAT 

CALL 'VECPRPC XHAT , BI GY > X 1K.»X2K,X 1KP1 , X2KP1 , TRMAT# THETA* ZETA# U# YM,YO 
+* TS.TEP> s XI I , P > NNUMER , N DEN 0 M j,N DAT A ) 

39 CONTINUE 

c' now ’Rave x -overbar, underbar and yo+ x-subii vector 

C S0LVE| FOR UPDATED P-VECT0K 

C LEAST 1 SQ. FIT FOR P-VECTOR 

C TRANSPOSE XHAT 

CALL MTRANS ( XHATT,XHAT , N DATA, NPRAM ) 

c cov=xhatt*xhat’ 

CALL ATMULC CO V,XHATT, XHAT * NPRAM* N DATA, NPRAM) 

C COVI= INVERSE (C0V ) 

CALL MATlNV(NPRAM,C0V,EPSl,EPS3,C0VI > 

C XHAT*BI GY REQUIRED TO GET P-VECTOR 

CALL Matmulcysol,xhatt,bigy,npram,ndata, JXK) 

C NOW READY FOR P-VECTOR '' 

C SAVE P, -VECTOR IN PBAC^ BEFORE GETTING NEW P-VECTOR 

DO 4.1 ;i=l, NDENOM 

1 1 = I +NNUMER+ 1 ' •' 

PBACKCil >=PC 1 1 ) 

41 CONTI NpE \ 

CALL MATMUL ( P, CO VI , YSOL,NPRAM, NPRAM, JXK ) 

c Compute? error .vector cerv), and sum sq. error cssqerj 

C NOW PREPARE TO COMPUTE MODEL OUTPUT ,'YM-USE NEW P-VEC 

CALL M^TEXP C TRMAT, P, TS TEP, NDENOM, NNUMER, NPRAM, D1 , D2, D3 ) * 

CALL . TRET AS t < THETA> P, TRMAT, NDENOM, NNUMER, NPRAM, D 1 , D2 > 

.CALL’ ZETAS < ZETA, P, THETA, TSTEP, NDENOM, NNUMER, NPRAM, D1 , D2 } 

C NOW,* USE VECPRP TO, GE\ XHAT 

C ALL VECPRP (XHAT ", BI GY, X 1 K, X2K, X 1KP 1 , X2KP 1 , TRMAT, THETA, ZETA, U, YM, © 
+ , TSTEP',‘Xl|l , P, NNUMER, N DENOM, NDATA ) 

C COMPUTE YM FOR POSSIBLE USE ON NEXT ITERATION 

DO .90' Kj= L, NDjAtA 
YM ( K ) = 0 • ' *' 

DO 100 1=1, NPRAM 
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YM ( K > = YM < K ) +XHAT ( K* I )*P( I ) 

100 CONTINUE 

DO 110 I=l*NDENOM 
1 1 =NNUMER+ 1 + 1 

YM ( K ) = YM C K ) - X HAT ( K* II >*PBACK(I ) 

110 CONTINUE 

90 CONTINUE 

YMCNDATA+ 1 )=YM(NDATA) 

C COMPUTE SUM S QUAKED ERROR 

SSQER=0. 

DO 120 K=1*NDATA 
EfiV ( K ) = Y 0 ( K ) - YM < K > 

120 S5QEK=$SQER+ERV ( K) *ERV ! K) 

C WRI TEC 7j 99 0) ( !P(I)* 1= 1*NPRAM ) » SSQER) 

NOINT=NINP- 1 

CALL WAVEHT C S WH> AT2* AT* A, P* TSTEP* NOINT* TRMAT, THET A* D2> ZETA* D1 * D3,P 
+ROB* NNUMER* NDENOM* NPRAM*NDATA) 

C WRITE ( 7* 505 > (SWK*AT2*AT* A) 

505 FORMAT C 4F 15*8) 

990 FORMAT (6F1 2 • 4 ) 

********.*TEST FOR STOPPING CONDITIONS ********** 

STOP WHEN SSQER BECOMES CONSTANT OR WHEN ITERN0>16 
FIRST CHECK FOR CONVERGENCE* THEN FOR ITERATION LIMIT 
ERCHK= SSQER- SSQERB 
ERCHK= DABS ( ERCHK ) 

SSQERD= SSQER 

IFCERCHK.LE. 1 .D-03) GO TO 200 
IF ( I TERNO . GT • 15) GO TO 250 
150 CONTINUE 

200 WHI TEC 7* 400) 

400 FORMAT!/ * K(S) PARAMETER VECTOR’) 

WRITE (7*401 ) (P(I )* I=1*NPRAM) 

401 FORMAT C6F1 2« 4 > 

WRITE! 7* 989) 

989 FORMAT! • PROB. DENSITY ’ > 

WRITE! 7* 99 0) (PROBd ) * I = 1 * NDATA ) 
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WRI 'I E< ( 7 , 968 ) 

988 FORMAT C \ SWH=, AT2=, AT=, A* : * > 

WRI TEC 7, 505') .CSUH, AT 2, AT,. A) 

WRITE<7:,402>' 

402 FORMAT < // ' RESIDUALS ’ ) 

• WRI TEC 7,401) CERVCI ), I=1,NDATA> 

• WRI TEC 7,403) SSQER 
GO TO 450 

403 F,ORMATC//' SSQER= * FI 2. 8 ) 

£50 WRIT EC T> 4 10) ‘ ' 

GO. TO 450 

420 WHI TEC 7,503) ' ‘ 

503 FORMAT C/// * CHECK INPUT DATA, H 1 3 COMPi NOT' ATTEMPTED') 
410 "FORMAT (.//'- *****dj VERGENCE***** ' ) ‘ ' ' ' “ 

450 CONTINUE 

RETURN 
' END 

* 
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SUBROUTINE MTRAMS C B, A, I RA, I CA> 

RETURNS TRANSPOSED) IN B 
IHA=#ROWS IN A 
I CA=#COLS. IN A 

. DIMENSION ACIRA, ICA),BCICA, IRA) 

DOUBLE PRECISION A,B 
DO 5 I = 1 , I RA 
DO 5 J=1,ICA 
BC J, I >=ACI, J> 

CONTINUE 

RETURN 

END 

-SUBROUTINE SKLMULC B, A, FACTOR, I RA, ICA) 
MULTIPLY A BY SCALAR=FACTOR, RETURN ■ IN B 
DIMENSION BCIRA, ICA), ACIRA, ICA) 

DOUBLE PRECISION B, A, FACTOR 
DO 5 1 = 1, IRA 
DO 5 J= 1,1 CA 
B( I , J) = FACTO h*AC I , J ) 

CONTINUE 

RETURN 

b;nd 


SUBROUTINE MATMULCC, A, B, I RA, I CRAB, I CB ) 

THIS ROUTINE COMPUTES A*B AND RETURNS IN C 
I HARROWS IN A 

ICRAB=0COLS. IN A,# ROWS IN B 
ICB=# COLS. IN B 

DIMENSION AC IRA, I CRAB ) , BC I CRAB, ICB),CCIRA, ICB) 
DOUBLE PRECISION A, B, C 
DO 5 1 = 1, IRA 
DO 5 J=1,ICB 
CC I, J) = 0. 

DO 5 K= 1,1 CRAB 
CCI, J)=CCI, J)+ACI,K)*B<K, J) 

5 CONTINUE 

RETURN 
END 
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SUBROUTINE ZETAS C Z ETA, P, THETA, TSTEP,NDENOM, NNUMER,NPRAM, C, AIN V) 
C THIS ROUTINE .COMPUTES ZETA=AINV*CTHETA-TSTEP*I ) 

DIMENSION THETACNDENOM,NDENOM>,PCNPRAM)>ZETA(NDENOM,NDENOM> 
DIMENSION. CCNDENOM,NDENOM ) , AINV CNDENQM, NDENOM ) 

DOUBLE PRECISION ZETA, P, THETA, TSTEP, C, AINV ' 

NDENO=NDENOM 
DO 5 I = 1, NDENOM 
DO 5 J= 1, NDENOM 
C C I ,_J) =THETAC I , J) 

IECI.EQ.J) CCI, J> = CCI, J>-T,STEP 
C SET UP AINU FROM P-VECTOR 

AINVC I , jW,0» 

IM= I - 1 

IFCIM.EQ.J) AINV C I , U) = 1 • 

NPICK=NNUMER+1 
NSEL=NPI CK+ J+ 1 
IFC J.EQ.NDENOM) GO TO 5 

IFCI .EQ. 1 ) AINVC I, J) = *-PCNSEL)/P<NPICK+1 ) 

5‘ CONTINUE 

AIN VC 1, NDENOM >='l./P : CNPICK+l) 

C NOW COMPUTE ZETA 

CALL' MATMUL CZET A, AINV, C,NDENO,NDENO,NDENO ) 

RETURN 

END’ 


SUBROUTINE MATADDCC, A,B, IRA, ICA) 

C ADD; A TO B, RETURN IN C 

D I MEN SI ON AC I RA, I CA) , B C I RA, I CA) , C C I RA, I CA ) 

DOUBLE. PRECISION C,A,B 

DO 5 1=1, IRA 

DO 5 J=1,ICA 

C.CI, J)=ACI, J5+BC-I, J) 

5 CONTINUE . 

RETURN ' 

END 



no. oono ono 


SUBROUTINE MATEXP < TRMAT, P, TSTEP, NDENOM, NNUMER,NPRAM, DM, FM, XM AT) 
THIS ROUTINE COMPUTES MATRIX EXPONENTIAL- CA+TSTEP) 

TSTEP=TIME STEP SIZE 
NDENOM=ORDER OF SQUARE MATRIX 
DIMENSION TRMAT (NDENOM, NDENOM), PCNPRAM) 

DIMENSION DM (NDENOM, N DENOM ) , FM ( N DENOM, N DENOM ) 

DIMENSION XMAT (NDENOM, NDENOM ) 

DOUBLE PRECISION I TRMAT, P, TSTEP, DM, FM, XMAT, TSPRME, DI VFAC 
INITIALIZE AS FOLLOWS: 

TRMAT TO IDENTITY MATRIX 
DM = TRMAT 

FM = FORMED FROM DENOMINATOR TERMS IN P-VECTOR 
TSPRME=TSTEP/256. 

DO 5 1=1, NDENOM 
DO 5 J=l, NDENOM 
IF(I.EQ.J) TRMAT ( I , J) = 1 • 

IFCI.NE.J) TRMAT Cl, <J> = 0* 

DMC I , J)=TRMAT C I , J) 

FM(I, J)=0. 

JJ=NNUMER+1+J 
IP1=I+1 

I F ( I . EQ. NDENOM) FM (I,J)=P< JJ) +TSPRME 
IF(IPl.EQ.J) FM(I , J>=TSPRME 
CONTINUE 
NDEN0=N DENOM 

COMPUTE POWERS OF ARGUMENT MATRIX FM AND ACCUMULATE INTO TRMAT 
RECAL THAT TRMAT IS INITIALIZED TO I -MATRIX 
DO 10 1=1,10 

CALL MATMUL (XMAT , DM, FM, NDENO , NDENO , NDENO ) 

DO 15 11=1, NDENOM 
DO 15 JJ=1, NDENOM 
DM (II, JJ)=XMAT <11, JJ) 

15 CONTINUE 

TEMP=FLOATCI ) 

DI VFAC= 1 . / DBLE C TEMP ) 

C MULTIPLY DM BY DIVFAC 
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CALL SKLMUL < DM* DM* DI VFAC*NDENO*NDENO ) 

C ACCUMULATE INTO THMAT 

■ CALL MATADD(TRMAT*TRMAT,DM*NDENO*NDENO) 

10 CONTINUE * 

C NOW RAISE TO' ROWER <TO CANCEL TSPRME EFFECT ABOVE) 

C -'USE DM ANDFM* FOR TEMP. STORAGE 

DO' 20 1 = 1*4 '■ 

CALL MAfMULCFM* THMAT* TRMAT*NDENO*NDENO*NDENO ) 

C NOW REPEAT TO GET TO CORRECT POWER 

g0 CALL MATMULC TRMAT* FM* FM*NDENO*NDENO*NDENO ) 

C TRANSITION MATRIX NOW STORED IN TRMAT 

• RETURN' 

END 


SUBROUTINE THETAS (THETA* P* TRMAT*NDENOM*NNUMER*NPRAM* C* AINV) 

C THI S ROUTINE COMPUTES THETA=INVC A) *(TRMAT-I ) 

DIMENSION THETA CNDENOM* NDENOM ) * P (NPRAM) * TRMAT (NDENOM* NDENOM ) 
DIMENSION C (NDENOM* NDENOM ) * AI NV (NDENOM* NDENOM ) 

DOUBLE PRECISION THETA* P* TRMAT* AINV* C 
NDENO=NDENOM 
DO 5 I = 1 * NDENOM 
DO 5 J=l* NDENOM 

C( I * J)=TRMAT( I * J) 

IF(I.EQ.J) C( I* J)=C(I* J)-l . 

C SET UP AINV FROM P- VECTOR 

AINJV(I* J> = 0.' 

IM-’I - 1 ' 

IFCIM.EQ.J) AI NV ( I * J> = 1 • 

NPI CK=NNUMER+ 1 
NSEL=NPICK+J+1 
IF(J.EQ. NDENOM) GO TO 5 

IF(I.EG.l) AINVU,* J> = -P(NSF,L)/P(NPICK+1 ) 

5 CONTINUE’ 

AI NV ( 1* NDENOM ) = 1 .'/P(NPI.CK+ I ) 

C NOW' COMPUTE THETA 

CALL MATMUL(THETA* AINV* C* NDENO*ND£NO* NDENO ) 

RETURN 

END 
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SUBKU U‘l INK VhCPKPCXHAT, B1 GY, X 1 K, X2K, X1W1 , X2KP 1 , TRMAT, THETA, ZETA,U 
+, YM, YO,TSTEP,XI X , P,NNUMER,NDENOM,NDATA) 

THIS ROUTINE PROPAGATES INITIAL STATE VECTOR THROUGH 
ENTIRE INPUT FORCING FUNCTION 

DIMENSION XHAT CNDATA, NDATA) , BI GY ( 1 6 ) , X1K C 8 ) , X2KC 8 ) , X 1KP 1 C 8 ) 
DIMENSION X2KP1 ( 8 ) , TRMAT (NDENOM, NDENOM ) , THETA (NDENOM, NDENOM) 
DIMENSION ZETACNDENOM, NDENOM), U( 1 7 ) , YM C 1 7 ) , YO C’l 7 ) , PC 12 ) 

DIMENSION XI I C 16 ) | 

DOUBLE PRECISION XHAT, El GY, X 1 K, X2K, X 1 HP 1 , X2KP 1 
DOUBLE PRECISION TRMAT, THETA, ZETA, U, YM, YO, P, XI I 
DOUBLE PRECISION TSTEP, DU, DY 
JXK=1 

INITIALIZE VECTORS FIRST 
DO 35 1=1, NDENOM 
XlKCI )=0. 

X1KPKI > = 0. 

X2KCI >=0. 

X2KP1 ( I )=0* 

CONTINUE 
DO 60 K=l, NDATA 

COMPUTE STATE VECTOR FOR H(S> SUB-I 
FIRST GET TRANSIENT COMPONET 

CALL MATMULCX1KP1, TRMAT, XI K, NDENOM, NDENOM, JXK) 

NO TE: X 1KP 1=X IK UPDATED ONE STEP BY STATE TRANSITION MATRIX 
NOW ADD FORCING TERM 
INUMP1=NNUMER+1 
DO 40 1=1, NDENOM 
61 X1KP1 (I )=X1KP1 ( I ) +THETA Cl, NDENOM >*UCK) 

DU=U(K+1)-U(K) 

DU=DU/TSTEP 

64 X1KP1CI )=X1KP1< I )+ZETA(I,NDENOM)*DU 

XHATCK, I ) =X 1 K C I ) 

C NOW SAVE X1KP1 FOR NEXT PASS THROUGH (K+ 1 ) 

XlKC I )=X1KP1 ( I ) 

40 CONTINUE 

C ***********NOW COMPUTE STATE FOR H(S) SUB-I I 

CALL MATMULCX2KP1, TRMAT, X2K, NDENOM, NDENOM, JXK) 
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XII <K)=0. ' 

BI GY ( K > = 0 • 

C XII < K ) IS ‘USED TO .TEMPO RAIRLY .HO^D (X-SUB-I I (TRANSPOSED) *PBACK) 

DO 50 IsIjNDENOM ' ‘ ' ‘ *’ ‘ 

X2KP1 < I >=X2KP1 < I ) +JHETAC I t NDENOM) *YM (K) 

DY=YM(K+1 )-YMCK) 

DY=DY/TSTEP ' ' 

X2KP1CI )=X2KP1 (I )+ZETA(I.,.NpENOM)*DY 
I2=NNUMER+I+1 ‘ : ’ ; " 

XHAT ( Kj 12) =X2K( I ) 

XII (K)=XII CK)+X2KCI )*PCI2) 

X2K(. I ,)=X2KP 1(1) 

50 CONTINUE 

BI GY ( K) = YO ( K) +XI I CK) '' 

60 CONTINUE 

RETURN 
END 
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SUBRO UTI ME WAV EH T C SWH, AT2, AT, A, P, TSTEP, NO INT, TRMAT, THETA, D2, ZETA.D 
+ 1, D3, PRO B, NNUMER, NDENOM, NPRAM, NDAT A) 

DIMENSION PCNPRAM),TRMAT(NDENOM, NDENOM), THETACNDENOM, NDENOM) 
DIMENSION PR0B(NDATA),D1 (NDENOM, NDENOM), D2 (NDENOM, NDENOM) 
DIMENSION D3(NDEN0M, NDENOM), ZETA(NDENOM, NDENOM) 

DOUBLE PRECISION P, TRMAT, THETA, D I , D2, D3, Z ETA 
DOUBLE PRECISION S WH, AT2, AT, A, TSTEP, BT, T, PARM 
THIS ROUTINE COMPUTES: 

SWH=H 1/3=S I GNI FI CANT WAVE HEIGHT 
AT2 = E(T**2),'E= EXPECTATION 
AT =E( T) 

A= I NTEGRAL ( H ( T ) -AREA UNDER PR0B. CURVE 
X=FLO ATCNO INT-4 ) 

BT=TSTEP*DBLE(X) 

T=BT/ 16. 

COMPUTE TRMAT 

CALL MATEXP C TRMAT, P, T, NDENOM, NNUMER, NPRAM, D1 , D2, D3 ) 

RAISE RESULT TO 16TH POWER 
DO 5 1=1,2 

CALL MATMULCD1, TRMAT, TRMAT, NDENOM, NDENOM, NDENOM) 

CALL MATMULCTRMAT, Dl, Dl, NDENOM, NDENOM, NDENOM) 

NOW COMPUTE THETA 

CALL THETAS ( THETA, P, TRMAT, NDENOM, NNUMER, NPRAM, Dl , D2 ) 

NOTE: AC **- 1 ) IS IN D2 FROM LAST CALL 
COMPUTE A: 

NNUMP 1 =NNUMER+ 1 
A=0. 

DO 10 1=1, NNUMP 1 
10 A=A+P( I )*THETA<I> NDENOM) 

C COMPUTE AT 

CALL SKLMUL ( Dl , TRMAT, BT, NDENOM, N DENOM ) 

CALL SKLMUL (D3, THETA, - 1 . , NDENOM, NDENOM ) 

CALL MATADDC Dl , Dl , D3, NDENOM, NDENOM ) 

C RECALL A**C-1) IS STORED IN D2 FROM THETAS CALL ABOVE 

CALL MATMUL<D3,D2,D1, NDENOM, NDENOM, NDENOM) 

AT=0. 
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DO 20 1 = 1* NNUMP1 
AT=AT+PCn*D3<I*NDEN0M> . 

COMPUTE AT2: 

PARM=BT*BT 

CALL SKLMUL < D 1 * TRM AT * PARM* N DENOM * NDENOM ) 

parm=-s.*bt' 

CALL SKLMUL ( D3* D2* PARM* NDENOM* NDENOM ) 

CALL MATMUL<ZETA* ! D3* TRM AT* NDENOM* NDENOM* NDENOM ) 

CALL MATADD’i D 1 * D 1 * Z ET A* N DENOM * N DENOM ) 

CALL MATMULCZETA* D2* THETA* NDENOM* NDENOM, NDENOM ) 

CALL SKLMUL C ZETA* ZET A* 2 • * NDENOM, NDENOM > 

CALL MATADDC D1 * ZETA* D.l* NDENOM* NDENOM) 

CALL MATMULCD3* D2* D'l * NDENOM* NDENOM* NDENOM ) 

AT2 = 0 * 

DO 30, 1 = 1* NNUMP 1 

AT2= AT2 + PC I ) *D3 C I * NDENOM) 

SWH='A^T2/A-(AT/A)*CAT/A) 

I F ( SWH • L T • 0.) GO TO 80 
SWH-DSQRT ( SWH) 

GO TO 9 0 
S WH= DABS ( S VH ) 

S WH= i- DS QRT < S l*iH ) 

SWH=6 • 25*5 WH/TSTEP 
SWH = DABSCSWH)*SVJH 

IFCS^HiLT. 0.) SWH=- DSQRT ( DABS ( SWH ) ) 

IFCSUH.GE. 0.) SWH=DSQRt.CSWH) . 

SWH=0.6*SWH 

COMPUTE PROB. DENS I TY ■ CURVE 
COMPUTE TRM AT FOR PROB. COMPUTATION 

CALL :MATEXP( TRMAT* P* T STEP* NDENOM* NNUMER*NPR AM* D1 * D2* D3 ) 
DO 40 1 = 1* NDENOM 
DO 40 J=l* NDENOM 
D1 ( I * J> = 0. 
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nid IFCI.EQ.J) D1 ( I * J)= 1 • 

DO 50 K= 1>NDATA 

CALL MATMULC D3* D1 > TRMATj NDENOMjNDENOM.»NDENOM> 
DO 60 I = 1 > NDENOM 
DO 60 J= 1 # NDENOM 
60 DKI, J> = D3d> J) ■ 

PRO B ( K) = 0 • 

DO 70 I*1,NNUMP1 
PARM=PCI >*DKI#NDENOM> 

7 tJ PROB<K)=PROBCK>+SNGL<PARM> 

PROB(K)=PROB<K)/SNGLCA) 

50 CONTINUE 

RETURN 
END 
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TT:*DK1 :>DMINV «FOR 

SUBROUTINE MAT I NV CNR, A 1 , EPS 1, EPS3, B 1 ) 

MATIN V COMPUTES THE INVERSE OF MATRIX A USING CROUT REDUCTION 
ANSWER RETURNED IN B 

E'PS INTOLERABLE ERROR OF RESIDUES IN SUBROUTINE CROUT 
EPS3=T0L • RELATIVE ERROR IN RESIDUES IN SUBROUTINE CROUT 
NR=NO • OF ROWS IN MATRIX A 
DIMENSION AC 13* 13>,B< 13, 13>,C( 13, 13), X( 12) 

DIMENSION A1CNR,NR),B1CNR,NR> 

DOUBLE PRECISION A 1 , EPS 1 , EPS 3, B1 , A, B, C, X 
DO I 1 = 1, NR 
DO 1 J=1,NR 

1 AC I , J) =A1 ( I , J) 

NSTO P=0 

DO 2 J=1,NR 
DC) 2 K=1,NR 

2 CCJ,K>=AC J,K> 

','NR1=NR+1 

DO! 6 K= 1 , NR 
DO' 4 J=. 1 , NR 

4 cc!j>nri')=0.' 
cc’k,nri>=i. , . 

CALL CROUT (NFL C, EPS 1, EPS3, X, NSTOP) 

: IF (NSTOP* EQ* 2 > GO TO 7 
DO '5 J=1,NR 

5 BC J,K)=X.( J,‘) 

6 CONTINUE, 

- GO.. TO 10 ' 

7 ' WRI TEC 7,50.6). < 

500 FORMAT (.{/// ’SINGULAR MATRIX ENCOUNTERED') 

10 CONTINUE’ ' ‘ * ‘ ' * ' 

DO 20 1 = 1, NR 
DO 20 J= 1 > NR 
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20 B1 Cl, J)=BCI, J) 

RETURN 
END 

SUBROUTINE CROUT CNR, B> EPS 1 , EPS3, T,NSTOP ) 

NR, Bj EPS 1 , EPS3 ARE NOT CHANGED BY THIS PROG* 

CROUT SOLVES FOLLOWING MATRIX EQ. : B CNR, NR > *BCNR+ 1 ) 

BCNRjNR) is to be inverted 

BCNR+ 1 > IS A KNOWN VECTOR 
SEE MATINV FOR DEFINITION OF EPS1,EPS3 
INVERTED MATRIX IS RETURNED IN VARIABLE T 
NSTOP IS FLAG FOR ZERO PIVOT ELEMENT 
DIMENSION A<12*12>*BC13#13>#X<12>jTU2)*R<12)jC(13> 
DOUBLE PRECISION A> Tj R>Cj EPS 1 , EPS3, Q, XX> RR 

DOUBLE PRECISION TEST 
501 FORMAT?/// ’MATRIX CANNOT BE INVERTED - SINGULAR’ > 
KP=0 
NSTOP=0 
NC=NR+1 
NGO=NC+l 

IFCNR.LE.3) NG0=6 
NN=NR- 1 
DO 2 1 = 1, NR 
RCI >=BCI,NC> 

2 T(I>=0 

3 DO 4 1 = 1, NR 
DO A J= 1 , NR 

4 AC I , J)=B( I , J) 

DO 5 N=1,NR 
ACN,NC )=RCN) 

C ( N ) = A C 1 , N ) 

5 CONTINUE 
CCNC )=A< 1,NC) 

DO 9 L= 1 , NN 
LO=L+ 1 

Q= DABS (AC 1, 1) )~DABS(ACLO, 1 ) ) 

IFCQ.GT. 0.) GO TO 9 
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I Sal •' . . ... 

IFCAC 1, 1 ).EQ. 0.) 60 TO 84 
DO 8 N=1>^C 
AC 1,N.)=.ACL0>N.) . 

A(LO.j>N)=;CCN) 

8 C<N^AUaN}. , ( 

9 CONTINUE'.. 

DO 11 I = 1 >.NR - 

ll.i AC I, I + 1,)=A( 1, I^1)/A(1j 1 > 

DO -50 .1= 1-,NN .• 

i s=i + t r >,’• 

-• NI=NR-.I •: 

DO i -3 0 K=1,NI 
KI=K+I. 

DO ; 3'0. L= !> I 

30 AC KLrl S ) = ACK,I * IS)-ACKI»L) *ACL* IS) 
DO 31 N=l,yNC 

31 CCN) = A,C IiSjN) 

LR=NK-IS 

IFCLR.LE.0) GO TO 39 
■ DO .38 L*= Tj.LR , 

LX-IS+L • • • 

. Q= DABS C A Cl S I S ) - DABS C A CLX* IS)) 
IFCQ..GE. 0.«) GO TO 38 
DO 35 N= 1 *!M0 ., 

AC IS*N ) -AC LX* N ) •. >. • 

A C LX j-N ) = 0 C N ) ■ 

35 CCN)=ACIS.,N) « .... . 

38 CONTINUE 

39 CONTINUE 

I F C AC I S> IS)* E0« 0 • ) , G 0. TO .84 
DO 45.J=-1j;NJ 
DO 40 M = 1 .» I 

I J=IS+0 

40 AC I S>I J) ~AC TS> I tJ[) ~AC I S» M ) *ACM> I J) 

45 A C I Sj> I J) =A CIS>IJ)/ACISj IS) 
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50. CONTINUE 

XCNR)=ACNR>NC> 

DO 60 KK=1jNN 

LL=NR-KK 

XX=0. 

1 1 =NR-LL 

DO 55 JJ=1,II 

MM=LL+JJ 

55 XX=XX+ACLL>MM>*X<MM> 

60 XCLL ) = AC LL.» NC ) -XX 

KTEST=0 
TEST=0. • 

DO 68 I=1>NR 
TCI ) = T C I )+X<I ) 

RR=0. 

DO 64 J=1,NR 

64 RR=RR+BCI, J)*XC J> 

RC I )=RC I ) -RR 

IFCDABSCXCI )'/T(I) >.GE.EPS3> KTEST= 1 
IFC A8SCRC I ) > • GE« EPS 1 ) TEST=I 
68 CONTINUE 

KP=KP+1 

IFC TEST .LE# 0. . AN D. KTEST . EQ . 0 -OR.KP. GT-NGO ) GO TO 86 
GO TO 3 
84 NST0P=2 

WRI TEC 7* 501 ) 

86 RETURN 

END 

* 
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INPUT NNUM* NDENOM* NINP-FO RMAT 
2 3 16 


313 


INPUT GATE SAMPLES 

! 


. * 

. 00 

• 

,.04 

• 

* 

05 


• 

.06 

• 

.08 



. 14 

.;25 

♦ 

32 


. 34 

. 50 



.68 

.82 

• 

85 


. 84 

.85 



.84 

!0 








NORMALIZED '.GATE 

: SAMPLES * 

* 





-0.0375 

0 .-0025 

0.0125 

0. 

0225 

0.0425 

0'. 1025 

0.2125 

0.2825 

0.3025 

0.4625- 

0. 6425 

0. 

7825 

0.8125 

0.8025 

0.8125 

0.8025 

0.8058 

0.8058 

0.8058 

.0.-8058 

.'0.8058 

0.8058 

0.8058 

0.8058 


0.8058 

HCS) parameter vector , 


60. 1 756 ’■ . 

-7.-8-903 

1. 1367 

-72.4544 

-36.7867 

- 

PROB. DENSI.T-y- * 
0.2029 

-0.-00,48 > 

0.3686 

0.9564 

1.4685 


1.7007 

■ 1,.;4 1.87 

0.9904 

0.5327 

0. 1432 

- 

-0.2396 

,-0.2372 

0.1539 

-0.0375 

0.0708 


0.1716 

■ 0. 1577 

0. 1155 

0.0616 

0.0112 

- 

-.-0.0438'. 
SWH=, AT2=>A.T=, 
9.06788551 ' 

A=,:- 

• . 0.36741901 

0. 

51712578 

0.84016753 



residuals 

-0.0375 

0.0534 

0.0225 

-0.0067 
J.0. 029,0- 
---04/020,5 

, -0.0118 
-0.0726 
-0.0179 

-0. 0179 
-0. 0438 
-0.0190 

-0.0221 

0.0146 

-0.0005 


• -0.0206 
-0.0283 

0> 0201 

0.0130 

0.0017 

-0*0105 



SSQER= 0.020^4843 
STOP -- 


5.7662 

1.7331 
0.1 195 
0. 1439 
0.0254 


0.0029 

0.0575 

0.0132 

0.0210 



Appendix B 

Plateau Region Sea State Behavior 


Theoretical analyses have shown that the plateau region, which contains 
stationary statistics, is characterized by a covariance function which is 
dependent on sea state [1]. A series of covariance calculations were made 
under contrasting sea state conditions to test this sensitivity. Figure B.l 
shows the results of these calculations; for the cases tested the change in 
covariance with sea state was found to be less than the uncertainty in the 
calculation. Based on this result, it- is concluded that sea state effects 
are analogous to a linear system convolution concept only in terms of a mean 
waveform effect and not in terms of the signal correlation properties. 
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